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SUMMARY 


Accurate  prediction  of  turbulence  inside  diesel  engine  combustion  chambers  is  one  of  the 
outstanding  obstacles  in  realistically  simulating  and  predicting  the  mixing  and  combustion  phenomena 
relevant  to  diesel  engines.  As  it  is  also  demonstrated  in  this  work  that  the  shortcomings  of  the  classical 
turbulence  models,  particularly  the  two-equation  models,  is  the  primary  factor  for  poor  predictions.  This 
study  has  been  performed  with  the  premise  that  at  least  some  of  these  turbulence  modeling  issues  can  be 
addressed  by  employing  the  large  eddy  simulation  (LES)  technique  to  predict  in-cylinder  flow  dynamics. 
In  the  LES  approach  the  most  important  (energetic),  large  turbulent  scales  are  resolved,  while  only  the 
effect  of  unresolved  small  scales  are  modeled.  For  the  present  simulations  a  combination  of  the  existing 
KIVA  family  (versions  2,  3,  and  3V)  of  codes  originated  from  Los  Alamos  National  Labs  and  a 
WVU/LES  code  developed  for  modeling  atmospheric  turbulence  were  utilized.  The  3-D,  unsteady  KJVA 
codes,  which  were  designed  specifically  for  internal  combustion  engines,  have  been  modified  for  LES  of 
turbulent  combustion.  The  WVU/LES  code  is  used  both  as  an  aid  in  understanding  the  effects  of 
swirling/rotation  on  turbulence  as  well  as  an  initial  benchmark  code  to  compare  how  well  KIVA/LES  is 
behaving  on  relatively  simple  flow  configurations. 

First,  to  identify  possible  problems  that  would  be  encountered  when  running  KIVA  in  LES  mode 
a  series  of  tests  were  performed  for  a  swirling  jet  flow.  This  flow  has  many  of  the  characteristics  of  the 
swirling  flow  induced  during  the  intake  stroke  to  aid  in  turbulent  mixing  and  combustion  inside  the 
engine  cylinder.  The  problem  was  simulated  using  KIVA-3  and  the  results  compared  favorably  to  those 
obtained  from  the  well-validated  WVU/LES  code.  Considering  the  differences  in  grid  resolution,  the  sub¬ 
grid  scale  (SGS)  model,  and  far-field  boundary  conditions,  reasonable  agreement  between  the  two  results 
was  observed.  Further,  a  new  SGS  model  applicable  to  strongly  rotational,  swirling  flows  was  developed 
using  the  WVU/LES  code,  which  has  significant  potential  for  IC  engine  modeling. 

Another  goal  of  this  research  was  to  improve  existing  RANS  (Reynolds  Averaged  Navier  Stokes) 
models  when  applied  to  IC-engines.  To  achieve  this  a  parallel  study  was  performed  on  a  benchmark 
flow,  which  resembles  that  of  a  motored  IC-engine  to  assess  the  performance  of  various  turbulence 
models.  As  expected  the  models  do  not  perform  uniformly  well  over  all  flow  regimes.  Significant 
differences  were  observed  among  various  models  as  the  engine  speed  increased.  A  new  hybrid  model  has 
been  proposed  which  tends  to  RANS  calculations  with  an  eddy  viscosity  model  and  to  LES  with  a 
Smagorinsky  SGS  model,  in  the  limit  of  coarse  and  fine  girds,  respectively.  This  is  especially  suitable  for 
engine  simulations  because  most  engine  simulations  are  inherently  three-dimensional  and  transient  as  in 
LES.  RANS  simulations  have  also  revealed  that  the  turbulence  length  scale  indicated  by  the  k-E  model 
provide  a  good  base  line  for  deciding  on  the  degree  of  resolution  needed  for  LES  simulations. 

As  a  first  step,  LES  of  the  compression  and  expansion  strokes  under  motored  engine  conditions 
were  considered.  Much  work  has  focused  on  the  compression  and  expansion  strokes,  due  to  the 
complications  introduced  by  the  valve  dynamics  during  the  intake  stroke.  A  comparison  of  the  results 
using  different  (SGS)  models  showed  that  the  simulations  with  a  SGS  model  produced  an  energy  cascade. 
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derived  from  a  spectral  analysis  of  the  fluctuating  data,  that  qualitatively  resembled  experimental  trends. 
To  achieve  this,  the  numerical  errors  needed  to  be  carefully  controlled  using  small  time  steps  with  second 
order  convection  schemes.  The  predicted  turbulence  intensity  was  low  but  showed  the  same  trend 
(roughly  a  linear  relationship  with  mean  piston  speed)  as  experiments.  This  lower  intensity  can  be  partly 
attributed  to  approximations  of  the  unknown  initial  conditions  inside  the  cylinder  which  may  be 
influenced  significantly  by  the  intake  stroke. 

During  the  compression  and  expansion  strokes,  turbulence  is  induced  primarily  by  the  geometry 
of  the  piston-cylinder  assembly,  in  particular  the  piston-bowl.  Most  of  the  turbulence  in  a  real  engine  is 
induced  during  the  intake  stroke.  For  the  next  level  of  realism,  simulations  of  the  intake  and  subsequent 
expansion  and  compression  strokes  were  performed  for  a  typical  two-valve  cylinder  assembly  using  the 
code  KIVA-3V.  A  multi  block  structured  mesh  was  set  up  in  agreement  with  the  experiments  of  Catania 
et  al.  (1996).  Two  sets  of  runs  were  performed  on  grids  with  different  resolution:  220,000  and  440,000 
nodes.  Here  too  a  study  was  performed  using  different  sub-grid  scale  and  RANS  turbulence  models 
including  k-e  model,  no  turbulence  model,  and  a  Smagorinsky  model.  The  study  showed  that  the  growth 
and  the  subsequent  decay  of  turbulence  during  the  intake  phase  predicted  with  the  Smagorinsky  subgrid- 
scale  model  agreed  well  with  experiments.  The  power  density  spectra  of  the  fluctuating  velocity 
components  were  compared  with  those  obtained  from  the  measurements.  This  comparison  showed  that  at 
least  some  of  the  inertial  range  dynamics  is  captured  in  the  simulations 

Throughout  the  studies  outlined  above  it  was  found  that  the  KIVA  codes  were  deficient  in  several 
areas  with  respect  to  performing  LES.  Improvements  to  both  the  efficiency  and  accuracy  of  the  KIVA 
code  were  implemented.  The  time  accuracy  of  the  code  was  made  fully  second-order  by  implementing  a 
combination  of  two-stage  Runge-Kutta  and  Adams-Bashforth  schemes  into  the  advection  phase.  Spatial 
accuracy  was  improved  by  introducing  an  advection  scheme  where  central  differencing  is  used  in  the 
momentum  equations  and  quasi-second  order  upwinding  is  used  in  the  scalar  equations.  Efficiency  of  the 
single  processor  version  of  KIVA  was  improved  by  up  to  20%  by  implementing  a  more  sophisticated 
preconditioning  scheme  in  the  pressure  solver.  The  overall  computational  performance  of  the  code  was 
improved  significantly  with  the  distributed-memory  implementation  of  KJVA-3  (KIVA-3/MPI)  based  on 
one-dimensional  domain  decomposition  for  parallel  execution.  This  current  version  is  capable  of  handling 
fixed  grids  on  any  given  number  of  processors.  KIVA-3/MPI  has  been  tested  on  several  hardware 
platforms  including  SGI  Origin  2000,  Cray  T3E,  and  Windows  NT  based  Intel  platforms.  The  code  also 
has  been  ported  to  the  20-node  Beowulf  type  DEC  Alpha  Linux  cluster  developed  for  this  project  at 
WVU.  With  K1VA-3/MPI  it  was  possible  to  attain  reasonable  turn-around  times  for  up  to  four  million 
grid  nodes  for  simple  flows  such  as  the  swirling  jet. 

With  the  groundwork  being  laid  for  performing  LES  of  motored  diesel  engine  cylinders,  a 
parallel  study  has  been  performed  focusing  on  development  of  combustion  models  appropriate  for  LES. 
The  first  such  study  involved  implementing  a  modified  form  of  an  Eddy  Breakup  (EBU)  turbulent 
combustion  model.  The  modification  allows  the  EBU  model  constant  to  be  dynamically  linked  to  the 
reaction  surface.  This  link  is  achieved  through  the  use  of  a  conditional  Probability  Density  Function 
(PDF)  which  is  evaluated  at  the  flame  sheet.  This  model  should  be  applicable  as  a  subgrid  scale 


iii 


turbulence/chemistry-interaction  model  for  LES  of  turbulent  combustion  in  diesel  engines.  The 
application  of  the  modified  EBU  model  to  a  D.  I.  diesel  engine  using  KJVA  III  has  shown  the  importance 
of  this  modification  to  the  modeling  of  chemistry/turbulence  interaction  in  engines.  Furthermore,  a 
validation  study  of  the  conserved  scalar/assumed  pdf  turbulent  mixing  model  has  been  performed.  This 
model  was  then  applied  to  an  experimentally  studied  diesel  engine  with  combustion.  The  computations 
indicate  that  this  approach  can  be  used  as  the  SGS  combustion  model  to  the  modeling  of  turbulent  mixing 
being  a  compromise  between  accuracy  and  computational  expense.  Given  the  importance  of  turbulence 
mixing  on  the  diffusion  phase  of  the  diesel  engine  combustion  process  an  additional  study  was  performed 
to  investigate  the  influence  of  the  finite-rate  chemistry  in  a  bluff-body  stabilized  diffusion  flame.  The 
computations  showed  that  both  constrained  equilibrium  and  flamelet  models  could  reproduce  the  macro 
features  of  the  flame.  However,  higher  order  correlations  were  under-predicted  at  the  jet  boundary.  Also, 
a  strained  laminar  flamelet  library  is  built  and  used  in  conjunction  with  CFD  simulation  of  a  turbulent 
H2/Air  diffusion  flame  to  investigate  the  finite  rate  chemistry  effects  on  the  thermal  NOx  formation.  The 
simulations  show  that  NO  is  a  species  that  is  highly  sensitive  to  the  scalar  dissipation  rate. 

Finally  spray  combustion  of  n-Tedradecane  is  modeled  using  a  7-reaction  Arrhenius  mechanism 
and  Magnussen  eddy-dissipation  combustion  model.  One  set  of  simulations  was  performed  for  the 
continuation  of  the  intake  flow  turbulence  with  valves  through  the  compression  stroke,  and  another  set 
was  for  the  compression  and  expansion  strokes.  These  three-dimensional  calculations  can  be  considered 
as  under-resolved  LES,  since  no  attempt  was  made  to  capture  thin  flamelets  via  a  rigorous  SGS 
combustion  model.  Nevertheless,  both  sets  of  calculations  showed  that  the  turbulence  intensity  increased 
significantly  with  the  onset  of  combustion.  Though,  there  is  peripheral  experimental  evidence  for  such  a 
phenomenon,  it  needs  to  be  investigated  further.  It  should  be  mentioned  that  extreme  difficulties  were 
encountered  during  spray  diffusion  computations.  With  LES  the  spray  penetration  seems  to  be  too  deep 
on  to  the  wall,  which  results  in  great  stability  problems  even  with  time  steps  as  low  as  1  .OE-09  seconds. 

On  another  front,  a  two-channel  Doppler  Global  Velocimetry  (DGV)  and  a  point  Doppler 
velocimetry  (PDV)  system  have  been  developed,  partially  funded  by  this  project,  and  have  been  used  to 
measure  the  mean  velocity  in  a  plane  for  a  rotating  wheel,  a  fully-developed  pipe  flow,  and  an  axi- 
symmetric  turbulent  jet  flow.  It  was  shown  that  mean  pipe  and  jet  DGV  data  compared  well  with  other 
experimental  data.  However,  RMS  velocity  fluctuations  did  not  compare  well.  With  minor 
improvements  both  systems  can  be  used  to  obtain  accurate  turbulence  data  for  validation  of  large  eddy 
simulations. 

As  a  result  of  this  study  a  total  of  37  papers  were  published  in  national  and  international 
conference  proceedings  and  journals.  The  project  involved  three  post  doctoral  fellows,  five  graduate 
students,  and  two  undergraduate  students.  In  addition,  the  research  team  has  actively  participated  in 
dissemination  of  information  by  attending  conferences  and  organizing  symposia  and  publishing  news 
letter  articles.  The  team  members  were  particularly  active  in  ASME  Internal  Engine  Combustion 
Division,  and  in  Multi-dimensional  Engine  Modeling  Sessions  organized  as  part  of  the  annual  SAE 
Congress  and  Exposition. 
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INTRODUCTION 


Turbulence  generated  during  the  intake  stroke  followed  by  the  turbulence  induced  by  the 
geometry  of  cylinder-piston  assembly  of  an  internal  combustion  (I.C.)  engine  during  the 
compression-expansion  stroke  greatly  influences  the  combustion  process  hence  fuel  efficiency 
and  pollutant  formation  (see  e.g.  Hey  wood,  1987).  The  classical  multidimensional  modeling 
approach  is  to  use  some  semi-empirical  turbulence  model  (e.g.  the  well-known  k-e  model)  to 
close  Reynolds/Favre  Averaged  Navier-Stokes  (RANS)  equations  of  motion.  These  models  have 
serious  shortcomings  in  that  they  introduce  crude  assumptions,  which  are  not  always  validated 
by  experiments.  One  such  assumption  is  the  local  isotropy  of  turbulence  inherent  to  most  eddy 
viscosity  models.  For  example,  in  the  standard  k-E  model  the  eddy  viscosity  is  isotropic,  and  the 
equation  for  the  turbulence  energy  dissipation  rate,  e,  is  modeled  using  the  local  isotropy 
assumption  (see  e.g.  Hossain  and  Rodi,  1982;  Speziale,  1998).  This  model  uses  only  one  length 
scale  and  a  time  scale  to  represent  a  whole  spectrum  of  turbulence  length  and  time  scales  ranging 
from  the  large  eddies  (proportional  to  the  size  of  the  combustion  chamber)  down  to  the 
Kolmogorov  scales.  Moreover,  most  of  these  models  require  that  the  turbulent  stress  tensor  be 
aligned  with  the  mean  rate-of-strain  tensor,  which  is  not  valid  for  many  complex  flows  such  as 
those  observed  in  IC-engines.  They  do  not  respond  to  rapid  rate  of  strain,  because  the  history 
effects  on  the  transport  of  Reynolds  stresses  are  not  accounted  for.  An  important  drawback  of  the 
k-e  model  is  that  it  gives  a  completely  unrealistic  representation  of  the  normal  stresses,  failing  to 
reproduce  the  strong  normal  stress  anisotropy  observed  essentially  in  all  shear  flows. 

Most  of  the  earlier  multi-dimensional  models  applied  to  I.C.  engines  have  used  RANS 
models  (see  e.g.  Amsden  et  al.,  1985;  Ramos,  1989;  Burgess  and  O’Rourke,  1993).  In  this 
approach  the  equations  governing  the  instantaneous  field  variables  are  averaged  over  a  time 
interval  which  is  relatively  short  so  as  not  to  render  total  elimination  of  the  time  variation  of 
important  unsteady  flow  events  but  still  long  to  filter  high  frequency  components  of  fluctuations. 
These  averaged  quantities  are  then  assumed  to  be  equivalent  to  ensemble-averaged  quantities  at 
the  same  point  in  space.  Extensive  testing  and  validation  studies  in  the  past  (Reitz  and  Rutland, 
1991;  Celik  et  al.,  1992;  Han  et  al.,  1996;  Smith  et  al.,  1997;  Bo  et  al.,  1997;  Yang  et  al.,  1998) 
have  shown  that  the  turbulence  models  introduced  semi-empirically  for  computation  of  higher 
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order  correlation  of  fluctuating  quantities  did  not  prove  universality,  and  the  success  of  these 
models  (not  to  be  underestimated)  is  case  specific.  In  Chapter  5,  as  an  example,  we  present  our 
RANS  calculations  using  various  turbulence  models  (Yavuz  and  Celik,  1999a,b)  with  some 
computations  performed  by  Yang  et  al.  (1998)  using  two  different  turbulence  models,  namely, 
the  standard  k-e  model  and  a  Reynolds  stress  model.  The  drastic  differences  between  the  results 
of  these  two  models  is  alarming;  particularly  in  temperature  distribution.  A  natural  question  is 
which  of  these  results  is  closer  to  reality?  We  believe  that  such  outstanding  turbulence  modeling 
issues  can  be  resolved  by  making  use  of  the  Large  Eddy  Simulation  (LES)  technique. 

LES  has  proven  to  be  a  very  promising  methodology  in  other  fields  of  applications  for 
prediction  of  turbulent  flow  quantities,  both  in  the  mean  and  instantaneously,  including  major 
statistical  properties  of  turbulence  (see  e.g.  (Reynolds,  1989;  Galperin  and  Orszag,  1993;  Ragab 
and  Piomelli,  1993,  Rodi  et  al.,  1997;  Piomelli,  1993  and  1998),  but  it  has  not  yet  found  its  due 
role  in  applications  to  in-cylinder  flows  for  I.C.  engines  (ICE).  Many  of  the  shortcomings  and 
assumptions  made  in  RANS  simulations  can  be  eliminated  by  employing  the  LES  technique 
which  has  been  widely  tested  for  simple  flows,  and  was  recently  extended  to  predict  complex 
turbulent  flows  encountered  in  realistic  engineering  applications.  This  technique  solves  the  three 
dimensional,  transient  Navier-Stokes  equations  after  applying  a  spatial  filtering  similar  to  time 
averaging.  Depending  on  the  filter  width,  which  is  usually  a  function  of  the  numerical  mesh  size, 
LES  can  capture  the  most  important  large-scale  fluctuations  in  the  flow  quantities,  leaving  only 
relatively  small  scales  to  be  modeled  empirically.  The  finer  the  grids  size,  the  better  is  the 
resolution  of  turbulence  scales,  and  the  less  important  are  the  modeling  assumptions.  At  the  sub¬ 
grid  level,  the  small-scale  structures  in  turbulence  are  believed  to  be  more  isotropic,  hence  the 
usual  eddy-viscosity  type  models  are  more  appropriate.  If  LES  can  be  successfully  applied  to 
engine  flows,  it  should  provide  insight  into  the  above  mentioned  controversies  related  to 
turbulence  models  and  enhance  our  understanding  of  in-cylinder  turbulence  generation  and 
control.  In  particular,  it  will  help  to  resolve  the  outstanding  questions  relevant  to  integral  length 
and  time  scales,  heat  and  mass  transfer  rates,  reaction  rates,  and  cycle-to-cycle  variations.  It  will 
also  provide  data  that  can  be  used  to  refine  the  classical  turbulence  models. 

As  mentioned  above,  in  the  last  decade  a  substantial  effort  has  been  put  into  advancing 
the  prediction  of  turbulence  in  non-reacting  flows  by  application  of  LES.  It  has  also  been  applied 
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with  considerable  success  to  turbulent  reacting  flows  (see  e.g.  Garrick,  1995;  Erlebacher  and 
Hussaini,  1993;  Menon  and  Jou,  1991;  Sykes  et  al.,  1990,  1994).  Reynolds  (1980)  suggested  that 
LES  was  probably  the  best  way  to  model  combustion  in  reciprocating  engines.  The  application 
of  LES  to  in  cylinder  turbulence  for  I.C.  engines  is  rare  because  of  complications  introduced  by 
compressibility,  complex  cylinder  and  valve  geometries  with  moving  boundaries,  and 
particularly  the  complication  due  to  turbulence-combustion  interaction.  One  of  the  early  attempts 
of  LES  for  engine  flows  was  presented  by  Naitoh  et  al.  (1992).  This  study  even  with  a  relatively 
coarse  grid,  and  first  order  Euler  time  marching  scheme  has  shown  the  great  potential  of  LES  for 
ICE  applications.  In  another  study,  Haworth  (1998)  reported  encouraging  results  using  a  finite 
volume  based  computer  code  for  predicting  the  ensemble  averaged  trends  for  an  experimental 
axisymmetric  engine  cylinder  motored  at  200  rpm.  A  comparative  study  by  Yavuz  and  Celik 
(1999)  indicates  that  classical  two-equation  turbulence  models  performed  well  for  this  case  in 
predicting  the  mean  flow  quantities  but  failed  to  predict  the  turbulence  quantities. 

This  report  presents  results  from  our  work  where  predictions  of  turbulent  fluctuations  and 
the  statistics  of  turbulence  quantities  inside  IC  engine  cylinders  were  attempted.  For  this  purpose, 
the  well  known  engine  simulation  code  KIVA  (a  widely  used  engine  simulation  code,  see  e.g. 
Amsden,  1993,1997)  is  used  with  a  special  precaution  to  keep  the  numerical  accuracy  at  a  higher 
level  than  usual  as  elucidated  in  section  2.3.  The  capabilities  of  the  KIVA  code  when  used  for 
LES  are  tested  against  benchmark  cases  such  as  lid-driven  cavity  flow,  swirling  and  non-swirling 
free  jet  flows.  The  code  is  then  applied  to  a  typical  engine  geometry  under  motored  conditions 
with  or  without  combustion.  The  predictions  are  qualitatively  compared  to  experiments  primarily 
to  investigate  trends  and  implications  of  certain  physical  mechanisms  in  turbulence  generation.  A 
thorough  quantitative  comparison  with  experimental  data  would  require  a  resolution  of  such 
issues  as  comparison  of  predicted  cycle-resolved  Favre  averaged  quantities  with  the  ensemble 
averaged  quantities  (which  are  usually  presented  in  experimental  works)  in  the  absence  of 
ergodicity.  The  ensemble  averages  contain  cycle-to-cycle  variations,  which  could  contribute  up 
to  30-40%  of  cycle-resolved  turbulence  (see  Valentino  et  al.,  1998;  Catania  et  al.,  1996,1997). 
But  simulations  are  possible  at  this  time  for  only  a  couple  of  cycles. 

In  the  course  of  this  study  it  was  necessary  to  investigate  the  numerical  accuracy  and  the 
efficiency  of  the  computer  code  used.  A  formally  second  order  time  marching  scheme  was 
devised  and  implemented  into  the  KIVA  code.  Preliminary  results  with  the  second  order  time 
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scheme  compared  well  with  the  first  order  scheme  provided  that  the  time  step  used  in  the  latter 
was  sufficiently  small.  Also,  a  parallel  version  of  the  same  code  was  developed  using  domain 
decomposition  with  MPI  (message  passing  interface)  for  non-reacting  flows  with  fixed 
boundaries.  The  scalability  of  the  parallel  version  has  been  shown. 

As  for  sub-grid  scale  (SGS)  models,  a  Smagorinsky  model,  and  a  one-equation  model 
were  used  for  the  flow  field.  For  combustion  the  well  known  eddy  break-up  model  was  used  as  a 
SGS  model  but  also  an  attempt  was  made  to  formulate  a  PDF  (Probability  Distribution  Function) 
based  combustion  model.  Preliminary  results  from  the  application  of  this  PDF  model  as  a  SGS 
combustion  model  are  presented. 
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Chapter  1 :  Computational  Procedure 


1 .1  About  the  Computer  Code 

The  calculations  were  performed  using  a  readily  available  computer  codes,  KIVA-3  and  - 
3V  (Amsden,  1993,1997;  see  also  Amsden  et  al.,  1985  and  1989)  with  modifications.  KIVA-3  is 
a  transient,  multiphase,  multidimensional,  arbitrary-mesh,  finite  volume  computational  fluid 
dynamics  (CFD)  program  widely  used  for  internal  combustion  engine  simulations.  It  is  based  on 
the  Arbitrary  Lagrangian  Eulerian  (ALE)  method  (Margolin,  1997).  The  numerical 
representation  of  the  convective/advective  terms  in  Eulerian  approach  lead  to  both  diffusion  and 
dispersion  errors,  which  create  difficulties  in  resolving  sharp  interfaces  in  flow  variables.  The 
ALE  method  seems  to  remedy  some  of  those  deficiencies.  The  concept  of  the  ALE  is  to  allow 
some  displacement  of  grid  nodes  during  one  time-step  interval.  A  more  detailed  discussion  on 
the  ALE  method  is  presented  in  Chapter  3.  The  method  is  typically  implemented  in  three  phases. 
The  first  phase  is  an  explicit  Lagrangian  update  of  the  equations  of  motion.  The  second  phase  is 
an  optional  implicit  phase  when  the  sound  waves  are  allowed  to  move  many  computational  cells 
per  time  step  if  the  material  velocities  are  smaller  than  the  speed  of  sound.  This  leads  to  a  greater 
computational  efficiency.  The  third  phase  is  the  remapping  (or  rezoning)  where  the  solution  from 
the  end  of  phase  two  is  mapped  back  onto  an  Eulerian  grid.  This  mapping  is  essentially  one  step 
of  a  conservative  advection  algorithm  (Margolin,  1997).  The  special  case  of  non-moving 
boundaries  where  the  cells  are  always  mapped  back  to  the  original  grid  is  the  Eulerian  limit.  The 
reader  is  referred  to  Amsden  et  al.,  (1993)  and  Hirt  et  al.,  (1997)  for  more  details  regarding  the 
ALE  method  and  the  KIVA  family  of  codes. 

1.2  Sub-grid  Scale  Models 

Favre  averaging  of  the  Navier-Stokes  equations  produces  turbulence  stresses  of  the  form 

Xij=-pu;uj  (1.1) 

Here  and  after  the  standard  tensor  notation  is  used  where  repeated  indices  denote 
summation. 
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A  similar  expression  is  obtained  when  the  governing  equations  are  filtered  in  the  large 
eddy  simulation  (LES)  approach  especially  when  a  box  filter  is  used.  The  Boussinesq  eddy 
viscosity  approximation  is  applied  to  determine  the  sub-grid  scale  (SGS)  stresses,  which  in  this 
study  are  expressed  as 

(^u)sGs  -|pk6ij  -|p,V-ii5j.  (1.2) 

where  is  the  mean  strain  rate  tensor,  and  6;^  is  the  Kronecker  delta  function.  The  extra  SGS 
stresses  arising  from  the  non-linearity  of  viscous  stresses  (see  Piomelli,  1998)  are  neglected. 

The  standard  k-e  model  that  is  implemented  in  KIVA-3  consists  of  the  following 
equations 


+  V  •  (pi^)  =  pkV  •  u  +  g :  Vu  +  V' 

dt  3  ~ 
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(1.4) 


with  the  model  constants  =  1 .44 ,  =  0.92 ,  =-1.0,  Pr;^  =1.0  and  Pr^  =1.3.  The  SGS 

heat  flux  vector  is  calculated  from  SGS  stresses  by  using  a  constant  turbulent  Prandtl  number, 
Prt  =0.74.  KIVA-3  includes  a  SGS  model,  which  utilizes  the  k-equation  given  above  along  with 
constraining  the  8  values  to  satisfy  the  following  inequality: 

LPr.(C..  -C.,) 

Lsgs  is  an  input  length  scale  whose  value  is  typically  some  measure  of  the  computational 
cell  dimension.  In  this  study  Lsgs  was  taken  to  be  equal  to  the  grid  related  length  scale,  Lo  = 

1  /O 

(AxAyAz)  ,  and  used  either  as  a  variable  length  scale  or  as  a  constant  with  a  typical  cell  size. 
This  model  will  be  referred  to  as  the  k-e  |  sgs  model. 

Modifications  were  made  to  the  KIVA  code  to  allow  a  classical  Smagorinsky  model 
where  the  eddy  viscosity  is  calculated  from 


^SGS 


(1.5) 
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V,  =(C.Lsos)“(s„sJ'‘ 


(1.6) 


The  model  constant  Cs  was  set  equal  to  0.2,  a  typical  value  used  in  the  literature  (see  e.g. 
Rodi  et  al.,  1997).  The  sub-grid  length-scale  was  estimated  as  an  average  computational  cell 
dimension  (see  Smith  et  al.,  1998). 

In  some  of  the  cases  presented  “no  SGS  turbulence  model”  was  used.  This  requires,  in 
most  instances,  some  numerical  diffusion  for  stability,  which  is  accomplished  via  the  use  of 
QSOU  (Quasi  Second  Order  Upwind)  scheme,  described  below,  instead  of  the  central 
differencing  (CD)  which  has  no  diffusion  error,  or  by  using  a  combination  of  CD  and  upwind 
differencing  heavily  biased  towards  CD.  It  should  be  noted  that  there  had  been  other  successful 
LES  studies  reported  in  the  literature  (see  e.g.  Kamamura  et  al.,^  and  Tamura  et  al.,  as  referenced 
in  the  review  by  Rodi  et  al.,  1997)  without  using  any  SGS  model.  As  such,  the  SGS  contribution 
to  the  eddy  viscosity  is  determined  by  the  numerical  diffusion  inherent  to  the  numerical  scheme 
used. 

The  above  turbulence  models  were  employed  in  conjunction  with  the  commonly  used 
law-of-the-wall  boundary  condition  which  is  implemented  in  KIVA-3.  The  basic  assumption 
here  is  that  the  interaction  between  the  modeled  near  wall  region  and  the  outer  region  is  weak  as 
observed  experimentally  by  Brooke  and  Hanratty  (1993)  for  a  range  of  flows. 

1.3  Numerical  Schemes 

Convective  terms  are  advanced  explicitly  in  time,  and  the  diffusion  terms  are  advanced 
explicitly,  implicitly,  or  semi-implicitly.  The  degree  to  which  the  diffusion  terms  are  implicitly 
discretized  is  based  on  a  combination  of  stability  and  efficiency  considerations.  Sub-time  steps 
(referred  to  as  sub-cycles)  are  taken  to  advance  the  convective  terms.  The  time  step  in  each  sub¬ 
cycle  is  based  on  Courant  stability  considerations.  The  convective  terms  are  advanced  with  time 
steps  that  are  the  same  or  smaller  than  what  is  used  for  the  diffusion  terms.  Though  the  overall 
time  accuracy  for  convection  terms  is  only  of  the  first-order,  the  global  time  step  is  based  on 
several  considerations  including  stability  and  several  accuracy  constraints.  To  increase  the  time 
accuracy  even  further  the  time  step  was  confined  to  the  interval  [1.5e-7  -  5.e-7  sec]  far  below  the 
Kolmogorov  time  scale,  which  is  in  the  order  of  lO"^  -10'^  seconds  for  a  typical  automotive  size 
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engine  ranging  at  speeds  1000-2000  RPM  (Celik  and  Yavuz,  1997).  This  precaution 
compensates  for  the  first  order  time  accuracy  of  convective  terms  and  guarantees  the  time 
resolution  required  for  LES.  For  these  reasons,  the  authors  strongly  believe  that  the  time 
accuracy  is  adequate  for  LES  despite  the  formal  order  of  accuracy,  as  it  is  demonstrated  below. 
The  dominating  errors  in  the  present  calculations  with  an  average  grid  resolution  of —1.0  mm  are 
due  to  spatial  discretization. 

All  spatial  derivatives  other  than  the  convective  terms  are  approximated  with  central 
differencing.  Spatial  accuracy  of  the  convective  terms  is  limited  to  several  choices  including 
first-order  upwinding  (FOUP),  not  used  in  this  study,  central  differencing  (CD),  a  user  defined 
scheme  that  combines  the  previous  two  choices,  and  a  flux-limiting  scheme  known  as  Quasi- 
Second-Order  Upwinding  (QSOU).  Both  the  first-order  upwinding  and  the  QSOU  schemes  are 
monotonic.  In  this  study  either  QSOU  or  CD  schemes  is  used.  More  details  about  the  numerical 
issues  can  be  found  in  Chapter  3. 


Further  justification  to  the  numerical  scheme  can  be  provided  by  an  order  of  magnitude 
analysis  of  the  truncation  error.  It  can  be  shown  that  in  our  case  the  normalized  truncation  error 
is  given  by 


T.E.  =  C 


VpAt 


3(aP,Ax^) 


where  a  typical  mean  piston  speed  is  Vp=  500  cm/sec,  dy  =  3  cm  is  the  valve  diameter  and  Ci,  C2 

are  assumed  to  be  well  behaving  functions  of  space  and  time  and  3{C|)s&(C2).  Then,  for 

At  =  1 .0x10'^  sec,  and  Ax  =  0. 1  cm. 

Term  I  VpdyAt 
Term  II  Ax^ 


which  shows  that  spatial  error  will  still  be  the  dominating  factor,  and  there  is  no  need  for  further 
refinement  of  the  time  scheme.  It  should  be  noted  that  there  are  other  successful  LES  studies  in 
the  literature  which  use  first  order  time  scheme  (see  the  work  of  Tamara  et  al.,  and  Onera  et  al.  as 
reported  by  Rodi  et  al.,  1997  in  an  LES  Workshop;  see  also  Naitoh  et  al.,  1992). 
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Chapter  2:  Rotational  Effects  in  an  LES  Model 

Swirl  is  often  an  important  part  of  the  flow  in  a  diesel  cylinder,  either  as  a  byproduct  of 
how  air  or  fuel  is  introduced  into  the  chamber,  or  by  special  arrangement  to  enhance  mixing  for 
combustion.  The  WVU/LES  code,  which  can  not  be  used  directly  to  simulate  flow  in  a  realistic 
cylinder  due  to  its  lack  of  sufficient  flexibility  in  boundary  conditions,  was  used  for  two  sub 
tasks  during  this  project.  One  was  to  simulate  swirling  jet  flow,  which  was  then  used  as  a 
benchmark  for  intercomparison  between  WVU/LES  and  KIVA.  The  second  was  to  derive  an 
LES  subgrid  model  capable  of  including  some  rotational  damping  effects.  This  later  task  was 
completed  in  coordination  with  two  other  grants  dealing  with  large  eddy  simulations  of  swirling 
flows  (a  NSF  grant  on  modeling  Tornado  dynamics  and  a  NASA  grant  on  modeling  aircraft 
trailing  vortices).  Results  of  these  two  tasks  are  described  briefly  in  this  chapter.  More  complete 
details  may  be  found  for  the  swirling  jet  in  Lewellen  et  al.  1998,  and  for  the  subgrid  rotational 
damping  in  the  appendix  to  Lewellen  et  al.  1999. 

In  an  LES  one  explicitly  simulates  turbulent  eddies  large  enough  to  be  resolved 
on  the  computational  grid,  and  models  the  effects  of  eddies  which  are  smaller.  If,  in  a  given 
problem,  one  can  resolve  the  scales  of  eddies  most  important  for  turbulent  transport,  then  the 
most  important  role  of  the  subgrid  model  is  simply  to  provide  an  energy  sink  so  that  energy  does 
not  pile  up  on  the  grid  scale.  In  such  cases  the  simulation  results  should  be  insensitive  to  the 
details  of  the  subgrid  model.  Sometimes,  however,  there  are  regions  in  a  simulation  where  the 
larger  scale  turbulence  is  suppressed  for  some  physical  reason,  and  the  subgrid  contribution 
becomes  relatively  more  important.  A  familiar  example  is  near  a  wall,  where  the  scale  of  the 
most  important  eddies  is  limited  by  the  distance  to  the  wall.  One  must  insure  in  such  cases  that 
the  subgrid  model  is  not  overly  dissipative,  i.e.,  that  any  physics  suppressing  the  resolved 
turbulence  also  suppresses  the  turbulence  represented  by  the  subgrid  model.  In  the  appendix  to 
Lewellen  et  al.  1999,  we  detail  how  the  most  important  effects  of  rotational  damping  were 
included  within  our  subgrid  model,  in  analogy  with  how  damping  in  a  stratified  flow  has  been 
treated  in  previous  LES  studies  of  the  planetary  boundary  layer  (e.g.,  Lewellen  et  al.  1996). 

The  subgrid  model  used  in  WVU/LES  is  based  on  a  subgrid  turbulent  kinetic  energy 
equation,  and  a  diagnostic  variable  A,  representative  of  the  characteristic  size  of  eddies 
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comprising  the  subgrid  scale  turbulence.  A  rotational  damping  term  was  added  to  the 
determination  of  the  subgrid  turbulence  length  scale.  A,  allowing  rotation  to  have  the  same 
limiting  influence  on  the  scale  as  that  represented  by  the  stratification  damping  term  in  previous 
atmospheric  versions  of  the  code.  The  new  algorithm  is: 

A=  Min[  0.65z,  c  Max[Ax,Ay,Az],  q/(2N),  q/(2^)]  (2.1) 

where  q  is  equal  to  the  square  root  of  twice  the  subgrid  turbulent  kinetic  energy,  N  is  the  Brunt  - 
Vaisala  frequency  determined  by  stratification,  and  ^  is  a  measure  of  the  rotational  radial 
oscillation  frequency.  In  a  pure  axisymmetric,  rotating  flow 

^2  =  2rr/r3  (2.2) 

where  T  is  the  product  of  the  tangential  velocity  times  the  radius,  r,  and  the  prime  represents  the 
derivative  with  respect  to  r.  An  energy  balance  argument  shows  that  §  ^  satisfies  the  analogous 

criterion  for  axisymmetric  disturbances  in  a  pure  axisymmetric,  rotating  flow,  as  satisfies  in 
stratified  planar  flow. 

Of  course,  in  a  general  3-D  flow  it  is  difficult  to  measure  the  local  value  of  TP/r^  about 
it's  local  center  of  rotation.  The  problem  was  to  find  a  function  to  approximate  this  in  the  general 
case  where  the  local  center  of  rotation  is  unknown  and  may  be  moving.  Such  a  desired  function 
in  Cartesian  coordinates  was  found  to  be: 

4  2  =  2(Vp-V,i )( Vp.V,i  -  Vp-VVi)/(  Vp-Vp)  (2.3) 

Where  V  is  the  velocity  vector  and  p  is  the  pressure.  In  the  pure  axisymmetric,  rotating 
flow  Eq.  2.3  reduces  exactly  to  Eq.  2.2.  We  note  in  passing  that  the  stabilizing  effects  of  strong 
radial  temperature  gradients  within  a  rotating  flow  can  be  included  in  the  same  level  of 

approximation  by  adding  the  term  (-  Vp*VT/To)  to  the  right  hand  side  of  Eq.  (2.3). 

The  implementation  of  this  term  in  the  LES  code  appears  to  work  well.  Its  primary 
influence  is  to  reduce  the  small-scale  turbulent  diffusion  in  any  coherent  vortex  cores.  This 
would  have  no  effect  as  long  as  one  has  the  luxury  of  resolution  much  finer  than  any  local  vortex 
core  size,  but  this  is  not  likely  to  be  the  case  in  any  diesel  chamber  simulation,  particularly  when 
wall  effects  are  dominant. 
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Results  of  high  resolution,  fully  3 -dimensional,  unsteady  simulations  (utilizing 
WVU/LES)  of  a  strongly  swirling  annular  jet  are  presented  in  Lewellen  et  al.  1998.  The  goal 
was  to  provide  high  enough  resolution  to  permit  the  dominant  turbulent  eddies  to  be  examined. 
Near  the  nozzle  the  jet  spread  is  approximately  axisymmetric  and  driven  by  small  scale  (relative 
to  the  jet  radius)  turbulence;  further  downstream  the  spread  rate  increases  as  the  turbulence 
becomes  characterized  by  a  small  number  (~  2  -  4)  of  localized  parcels  of  concentrated  angular 
momentum  spiraling  radially  outwards.  Mean  velocities,  turbulent  variances,  and  joint 
probability  distributions  are  compared  with  laboratory  experimental  data  previously  obtained  by 
Holzapfel  (1996).  The  simulated  jet  dynamics  is  found  to  be  quite  sensitive  to  low- wave  number 
perturbations  of  the  velocity  distributions  issuing  from  the  nozzle  with  differences  between 
laboratory  and  numerical  data  easily  masked  by  the  uncertainty  in  this  boundary  condition. 
Definitive  comparisons  between  simulations  and  laboratory  results  would  require  quite  detailed 
specification  of  the  fluctuations  within  the  flow  from  the  nozzle.  This  inherent  sensitivity  of  the 
swirling  jet  to  details  of  the  inlet  conditions  suggests  that  it  should  not  be  used  as  a  standard  for 
fine-tuning  closure  model  parameterizations.  However,  qualitative  conclusions  can  be  drawn  by 
comparison  of  these  results  with  the  KIVA  predictions  (see  Chapter  3). 

Figure  2.1  shows  an  instantaneous  view  of  the  jet  vertical  velocity  in  a  plane  cutting 
through  the  center  of  the  annular  swirling  jet  utilizing  1.5  x  10^  grid  points  distributed  over  a  grid 
stretched  in  all  three  directions.  The  axisymmetric  average  over  a  few  flow  through  times  is 
superimposed  as  contours  over  the  instantaneous  view.  It  is  these  averages  which  are  compared 
to  the  laboratory  data.  Comparison  tests  (Chapter  3)  with  KIVA  utilize  both  these  averages  and 
instantaneous  distributions  for  a  qualitative  comparison  of  the  turbulent  fluctuations. 
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Figure  2.1  Vertical  velocity  distribution  obtained  from  WVU/LES  simulation  of  an  annular 
swirling  jet.  Time  averaged  contours  at  intervals  of  5  m/sec  are  superimposed  over  an 
instantaneous  view  represented  by  shading. 


Chapter  3:  Validation  of  KIVA  and  the  ALE  Method  for  LES  Application 


3.1  Background 

There  are  several  popular  solution  strategies  that  are  employed  in  performing  LES.  The 
most  effective  appear  to  be  spectral  methods.  This  is  due  to  the  high  accuracy  and  energy 
conserving  properties  associated  with  these  methods.  However,  spectral  methods  are  not  as 
flexible  as  traditional  finite  volume  methods  due  to  difficulties  in  applying  them  to  non-uniform 
grids,  a  necessity  for  representing  complex  geometries.  Finite  volume  methods  are  easily 
applied  to  non-uniform  grids,  and  as  a  result  can  represent  complex  geometries.  The  most  severe 
drawback  to  using  finite  volume  methods  to  perform  LES  is  related  to  accuracy,  a  problem  made 
worse  when  non-uniform  grids  are  used.  When  finite  volume  based  codes  are  applied  in  the  LES 
approach  accurate  non-diffusive  convective  differencing  schemes  are  required.  Fine  grid 
resolution  is  also  necessary.  The  later  problem  requires  that  a  relatively  more  efficient 
implementation  has  to  be  used.  The  ALE  (Arbitrary  Lagrangian-Eulerian)  method  may  be 
computationally  more,  efficient  than  other  finite  volume  implementations,  especially  those  based 
on  variants  of  the  SIMPLE  algorithm  (Patankar,  1980).  In  addition,  there  is  a  possibility  that  the 
convective  modeling  procedure  utilized  in  the  ALE  method  may  be  more  accurate  than 
traditional  approaches.  This  will  be  discussed  in  more  detail  in  the  next  section. 

Swirling  flows  are  of  great  interest  to  the  engineering  community  due  to  their  widespread 
application  in  engineering  systems.  They  are  of  particular  interest  to  both  the  gas  turbine  and  the 
internal  combustion  engine  conununities.  In  the  later  example  significant  vorticity  is  induced 
during  the  intake  stroke  to  aid  in  turbulent  mixing  and  combustion  inside  the  engine  cylinder. 
One  of  the  keys  to  the  successful  simulation  of  internal  combustion  engine  flows  is  the  accurate 
modeling  of  the  effects  of  rotation  on  the  mean  flow  field.  However,  the  simulation  of  the  intake 
stroke  of  an  IC  engine  is  too  complicated  of  a  problem  to  use  as  a  test  case.  It  is  desired  to 
simulate  a  flow  that  has  similar  swirling  properties  to  that  of  flows  commonly  encountered  in 
engineering  systems,  yet  is  simple  enough  such  that  detailed  experimental  data  exists  for 
comparison.  A  good  choice  is  the  case  of  an  annular  swirling  jet  such  as  that  recently  studied  by 
Holzapfel  (1996). 
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The  main  objective  of  this  study  is  to  assess  the  LES  capabilities  of  the  ALE  method.  In 
addition,  the  same  assessment  is  desired  of  the  commonly  used  computer  code  KIVA-3  (Amsden 
et  al,  1989,1993)  which  is  based  on  the  ALE  method.  This  assessment  will  be  performed  by 
simulating  a  free  three-dimensional  swirling  jet.  This  problem  is  complicated  mainly  due  to  the 
presence  of  rotation.  It  is  well  known  that  many  traditional  turbulence  models  (e.g.  standard  k-s 
model)  fail  in  the  presence  of  rotation  due  to  the  isotropic  Boussinesq  eddy  viscosity 
approximation.  Many  add-hoc  modifications  have  been  made  to  the  k-e  equations,  specifically 
to  the  8  equation,  to  capture  the  effects  of  rotation  in  the  traditional  Reynolds  Averaged  Navier- 
Stokes  (RANS)  models  (Bardina  et  al,  1985;  Launder  et  al,  1977;  Sloan  et  al,  1986).  However, 
when  similar  models  are  applied  as  sub-grid  scale  (SGS)  models  in  conjunction  with  Large  Eddy 
Simulation  (LES),  there  is  some  discrepancy  in  the  literature  as  to  whether  there  is  a  need  for 
such  modifications.  Some  researchers  believe  that  extra  rates  of  strain  such  as  rotation, 
streamline  curvature,  and  stratification  affect  the  large  scales  much  more  than  the  small  scales 
(Ferziger,  1993).  If  this  is  true,  then  it  should  be  possible  to  use  SGS  models  without  major 
modifications  for  flows  where  these  extra  rates  of  strain  are  significant.  Previous  investigations 
(Smith  et  al.,  1997,  1998)  have  shown  promise  in  capturing  a  swirling  flow  field  with  relatively 
simple  SGS  turbulence  models  when  sufficient  grid  resolution  is  used. 

In  this  study,  a  simple  Smagorinsky  SGS  model  is  used  in  addition  to  relatively  fine  grid 
resolution.  The  flow  field  obtained  from  the  simulations  is  compared  to  experimental  data  for 
the  same  flow  (Holzapfel,  F.,  1996).  The  turbulence  quantities  are  compared  to  assess  the  LES 
capabilities  of  the  ALE  method.  Previous  work  on  this  particular  flow  situation  has  been 
performed  using  a  widely  tested  LES  code  (Lewellen  et  al.,  1996,  1998),  hereafter  referred  to  as 
the  WVU-LES  code.  The  results  and  experience  gained  from  applying  the  WVU-LES  code  to 
this  particular  problem  were  used  as  a  guide  for  this  work.  In  addition  to  comparing  to  the 
experimental  data,  the  results  obtained  in  this  work  are  compared  to  those  obtained  with  the 
WVU-LES  code.  Qualitative  duplication  of  the  WVU-LES  results  using  KIVA  is  desired  and 
would  be  further  validation  of  the  KIVA  code  for  use  in  LES. 

The  swirling  jet  test  case  will  also  be  used  as  validation  of  a  newly  implemented  second- 
order  Runge-Kutta  time  integration  scheme. 
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3.2  Computational  Issues 


The  KIVA-3  code  is  based  on  the  Arbitrary  Lagrangian  Eulerian  (ALE)  method.  The 
ALE  algorithm  combines  features  of  Lagrangian  and  Eulerian  representations.  In  a  Lagrangian 
simulation,  the  mesh  moves  with  the  local  velocity  of  the  fluid.  An  advantage  of  the  Lagrangian 
representation  is  that  in  this  reference  frame  convective/advective  terms  vanish.  However, 
Lagrangian  meshes  tend  to  tangle  and  in  general  can  not  represent  large  deformations  that  result 
from  shear  and  vorticity.  In  an  Eulerian  simulation  the  mesh  is  fixed  in  space  and  fluid  moves 
from  one  cell  to  another.  The  Eulerian  methodology  does  not  lead  to  tangled  meshes.  However, 
the  numerical  representation  of  the  convective/advective  terms  leads  to  both  diffusion  and 
dispersion  errors  that  create  difficulties  in  resolving  sharp  interfaces.  It  is  these  errors,  especially 
numerical  diffusion,  which  make  finite  difference  applications  in  LES  difficult. 

The  Lagrangian  and  Eulerian  methodologies  are  two  special  cases  of  mesh  motion.  The 
concept  of  the  ALE  method  is  that  mesh  motion  can  be  chosen  arbitrarily.  The  method  is 
typically  implemented  in  three  phases.  The  first  phase  is  an  explicit  Lagrangian  update  of  the 
equations  of  motion.  The  second  phase  is  an  optional  implicit  phase  that  allows  sound  waves  to 
move  many  computational  cells  per  time  step  if  the  material  velocities  are  smaller  than  the  fluid 
sound  speed.  This  allows  for  greater  computational  efficiency.  The  third  and  final  phase  is  the 
remap  (or  rezone)  phase  where  the  solution  from  the  end  of  phase  two  is  mapped  back  onto  an 
Eulerian  grid.  This  mapping  is  essentially  one  step  of  a  conservative  advection  algorithm 
(Margolin,  1997).  The  special  case  of  no  boundary  motion  where  the  cells  are  always  mapped 
back  to  the  original  grid  is  the  Eulerian  limit.  The  reader  is  referred  to  Amsden  et  al.  (1989, 
1993)  and  Hirt  et  al.  (1997)  for  more  details  regarding  the  ALE  method  and  the  KIVA  family  of 
codes. 

The  underlying  computational  methodology  of  the  ALE  method  differs  considerably 
from  more  popular  methods.  However,  from  the  user’s  point  of  view  the  ALE  methodology  is 
not  that  different.  Time  advancement  is  similar  to  many  codes  that  utilize  the  LES  methodology 
where  the  convective  terms  are  advanced  explicitly  and  the  diffusion  terms  are  advanced 
explicitly,  implicitly,  or  semi-implicitly.  The  degree  to  which  the  diffusion  terms  are  implicitly 
discretized  is  based  on  a  combination  of  stability  and  efficiency  considerations.  It  should  be 
noted  that  the  user  can  set  the  implicitness  factor.  Sub-time  steps  (referred  to  as  sub-cycles)  are 
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taken  to  advance  the  convective  terms.  The  size  of  the  time  step  for  each  sub-cycle  is  based  on 
stability  considerations.  Although  the  overall  time  accuracy  is  only  first-order  accurate,  the 
global  time  step  is  based  on  several  considerations  including  stability  and  several  accuracy 
constraints.  The  convective  terms  are  advanced  with  time  steps  that  are  the  same  or  smaller  than 
what  is  used  for  the  diffusion  terms.  Due  to  these  reasons,  the  present  authors  believe  that  the 
time  accuracy  is  quite  good  despite  the  formal  order  of  accuracy.  To  further  insure  time 
accuracy,  the  diffusion  terms  were  advanced  in  a  second-order  Crank-Nicolson  approach  with  an 
implicitness  factor  of  0.5.  The  minimum  number  of  subcycles  taken  for  the  convective  terms 
was  set  at  10.  This  reduces  the  leading  temporal  error  by  an  order  of  magnitude.  It  should  be 
mentioned  that  second-order  time  accuracy  is  still  not  achieved.  However,  the  implementation 
outlined  above  should  yield  time  accuracy  that  is  adequate  for  this  problem  given  the  accuracy 
restraints  on  the  time  step.  Also,  the  size  of  the  time  step  that  was  used  for  each  sub-cycle  was 
much  smaller  than  dictated  by  the  accuracy  constraints.  In  addition,  a  fully  second-order,  time 
accurate  Runge-Kutta  time  marching  scheme  is  implemented  and  tested.  It  should  be  noted  that 
even  with  second-order  accuracy,  small  time  steps  need  to  be  taken  to  resolve  the  turbulence 
spectra  (Choi  and  Moin,  1994). 

The  QSOU  scheme  used  in  the  KIVA  code  approaches  second-order  accuracy  in  the 
absence  of  steep  gradients.  However,  when  it  is  desired  to  resolve  turbulence  in  combination 
with  subgrid-scale  models  CD  is  typically  used  (Piomelli,  1998;  Mittal  and  Moin,  1997)  to 
prevent  numerical  diffusion  from  entering  the  calculation.  Special  care  must  be  taken  in  applying 
CD  in  a  compressible  code.  Small  dispersion  errors  can  be  catastrophic  to  a  compressible  code 
such  as  KIVA.  For  a  LES  application  the  added  turbulent  eddy  viscosity  is  usually  enough  to 
ensure  that  the  effective  Peclet  numbers  are  small.  Thus  the  dispersive  errors  are  small  in  the 
momentiun  equations.  However,  even  small  dispersion  errors  lead  to  problems  when  the 
advective  schemes  are  applied  to  scalar  transport  equations.  For  example,  the  internal  energy 
equation  is  particularly  sensitive  to  dispersion  errors  and  can  lead  to  non-physical  density  and 
temperature  values.  These  quantities  are  always  solved  for  in  the  presently  used  ALE  approach 
regardless  of  whether  the  code  is  applied  to  a  compressible  or  incompressible  situation.  Small 
oscillations  in  these  variables  can  grow  and  lead  to  divergence  of  the  solution.  For  the  problem 
studied  here  variation  in  temperature  and  density  should  be  small  and  have  little  or  no  effect  on 
the  flow  field.  To  ensure  that  this  is  the  case  a  modification  was  made  to  the  code.  Whenever 
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CD  is  used  for  the  convective  terms  in  the  momentum  equation,  QSOU  is  used  for  all  scalar 
variables.  This  appears  to  have  negligible  effect  on  the  simulations  since  both  temperature  and 
density  remain  relatively  constant  (less  than  2%  variation).  This  is  an  excellent  result  for  what  is 
normally  a  compressible  code. 

By  combining  the  Lagrangian  and  Eulerian  approaches  in  the  ALE  method  it  is  hoped 
that  some  of  the  discretization  errors  associated  with  the  approximation  of  the  convective  terms 
will  be  eliminated.  This  would  make  the  ALE  approach  attractive  as  a  LES  methodology  since 
the  major  problem  facing  standard  finite  difference  approaches  is  the  error  resulting  fi’om  the 
convective  discretization. 

For  the  present  test  case  a  standard  Smagorinsky  model  was  implemented  as  described  in 
Section  1.2. 

Cs  was  taken  as  0.2,  a  typical  value  used  in  the  literature  for  isotropic  turbulence.  Lsos  is 
an  input  length  scale  whose  value  is  typically  some  measure  of  the  computational  cell  dimension. 
Here  it  was  calculated  from 

Lsos  =  (^xAyAz)^  (3.1) 

3.3  Application 

The  geometry  of  the  free  annular  jet  was  defined  to  be  the  same  as  that  experimentally 
studied  by  Holzapfel  (1996).  Laboratory  conditions  were  duplicated  as  close  as  possible  in 
formulating  the  numerical  simulation.  The  computational  domain  was  defined  as  a  three- 
dimensional  prism,  with  the  jet  inlet  at  the  bottom,  and  the  outlet  at  the  top.  The  jet  inlet  was 
defined  by  an  annulus  that  ranged  from  R=2.7-5.3  cm  on  the  xy-plane  at  the  bottom  of  the 
domain.  The  width  of  the  domain  was  defined  at  100  cm.  The  axial  length  was  set  at  100  cm, 
which  translates  into  40  aimulus  widths.  The  far  field  boundary  on  the  left,  right,  front  and 
derriere  sides  of  the  domain  consisted  of  setting  the  radial  entrainment  velocity.  This  was 
obtained  from  extrapolating  the  values  from  the  experimental  data  while  at  the  same  time 
satisfying  global  continuity.  It  is  believed  that  the  domain  width  is  of  sufficient  size  such  that 
any  errors  in  the  far-field  boundary  condition  will  have  little  effect  on  the  jet.  The  axial,  radial, 
and  swirl  velocities  at  the  inlet  are  given  by  Holzapfel  (1996).  The  same  profiles  were  measured 
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at  different  axial  locations  through  the  use  of  a  hot-wire  technique.  The  simulation  was 
performed  by  assigning  the  velocity  profiles  at  x=0.0,  and  then  letting  the  jet  develop  as  it 
propagated  downstream.  The  inlet  plane  itself  consisted  of  a  flat  plate,  the  sole  purpose  of  which 
was  to  guarantee  radial  entrainment  near  the  inlet.  The  log-law  boundary  layer  approximation 
was  used  for  this  boundary.  This  approximation  should  be  adequate  given  the  fact  that  the 
turbulence  in  this  flow  field  is  generated  at  the  inlet.  The  Reynolds  number  for  the  jet  based  on 
the  average  axial  velocity  at  the  inlet  and  the  width  of  the  jet  annulus  is  approximately  20,000. 

The  radial  plane  was  discretized  using  a  78x78  non-uniform  grid  as  is  shown  in  Figure 
3.1.  The  axial  (z)  direction  was  discretized  using  82  vertices.  The  computational  grid  thus 
formed  consisted  of  approximately  500,000  nodes.  Since  the  computational  domain  was  defined 
as  a  prism,  a  Cartesian  grid  system  could  be  used.  This  allowed  more  control  over  the  numerical 
errors  due  to  the  discretization.  Though  the  choice  of  a  Cartesian-type  grid  helps  to  minimize 
discretization  errors,  it  introduces  errors  in  approximating  circular  shapes,  in  this  case  the 
aimular  inlet.  To  better  resolve  the  aimular  jet  at  the  inlet,  the  grid  was  contracted  near  the  center 
of  the  xy-plane.  The  contraction  was  kept  at  slightly  less  than  3%.  It  has  been  found  that  when 
CD  is  used  grid  stretching  in  excess  of  3%  introduces  grid  to  grid  de-coupling  that  is  manifested 
by  oscillations  (Mittal  and  Moin,  1997).  This  constraint  was  relaxed  in  the  far  field  since  this 
region  should  have  little  effect  on  the  flow  field.  Even  with  grid  clustering  the  annulus  was  not 
perfectly  defined.  However,  this  error  helps  in  introducing  small-scale  perturbations  in  the  flow 
field.  These  perturbations  are  a  result  of  the  combination  of  strong  shear  zones  that  exist  near  the 
inlet  and  the  asymmetries  that  exist  in  the  grid.  This  helps  initialize  the  turbulent  jet. 
Introducing  random  perturbations  (10%)  to  all  three  components  of  the  inlet  velocity  further 
initializes  the  turbulent  field.  It  should  be  noted  that  the  turbulence  intensities  reported  by  the 
experiments  were  much  higher  than  this. 

The  axial  direction  was  also  defined  with  a  non-uniform  grid  (Figure  3.2)  and  consisted 
of  82  vertices.  The  vertices  were  placed  uniformly  for  the  first  10-cm  of  the  jet  at  an  interval  of 
0.5-cm.  This  allowed  for  good  resolution  for  the  first  jet  diameter  where  the  jet  is  particularly 
sensitive  to  the  inlet  conditions.  The  grid  was  allowed  to  expand  at  1.7%  from  z  =  10  cm  to  z  = 
20  cm.  The  grid  was  expanded  at  a  rate  of  2.7%  from  z  =  20  cm  to  z  =  52  cm.  The  region  of 
interest  is  confined  to  z<  40-cm  since  this  is  where  experimental  data  is  available.  With  this  in 
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mind  the  expansion  was  increased  to  10%  for  the  remaining  48-cm.  It  is  this  last  region  where 
the  effects  of  the  outflow  boundary  conditions  are  expected  to  be  nullified.  It  should  be  noted 
that  the  Smagorinsky  model  introduces  enough  eddy  viscosity  such  that  the  effective  cell  Peclet 
numbers  are  small.  This  helps  prevent  the  grid-to-grid  decoupling  that  occurs  when  CD  is  used. 

Two  different  outflow  boundary  conditions  are  available  in  KIVA-3.  The  first  is  a  zero 
slope  condition  for  the  velocity  gradients.  The  second  involves  specifying  the  pressure  at  the 
outflow.  Since  this  flow  is  highly  convectively  dominated,  the  outflow  boundary  condition 
should  not  have  significant  effect  on  the  upstream  velocity.  In  fact  previous  work  on  the  same 
problem  using  the  WVU-LES  code  (Lewellen  et  al.,  1998)  has  confirmed  this.  Also,  it  is 
believed  that  the  domain  is  long  enough  so  as  to  eliminate  the  influence  of  the  outflow  boundary 
condition  on  the  flow  region  of  interest.  The  zero  slope  boundary  condition  was  used  as  is 
commonly  employed  for  this  type  of  flow  (Lewellen  et  al.  1998;  Xia  et  al.,  1997). 

The  swirl  strength  is  characterized  by  the  jet  inlet  swirl  number  So .  So  is  defined  as  the 
ratio  of  angular  to  axial  momentum  flux  in  the  nozzle  divided  by  the  outside  nozzle  radius.  In 
the  experiments  an  approximate  swirl  number  was  calculated  based  on  the  swirl  generator 
geometry.  The  free  jet  case  that  was  simulated  in  this  work  had  a  swirl  ratio  (So)  equal  to  0.95  as 
reported  by  Holzapfel  (1996). 

3.4  Results 

3.4.1  Numerical  Data  Acquisition 

The  jet  was  allowed  to  develop  from  an  initial  state  of  zero  velocity.  This  consisted  of 
simulating  the  jet  for  approximately  five  flow-through  times  using  relatively  large  time  steps.  A 
flow-through  time  is  defined  as  a  time  scale  representing  the  average  length  of  time  required  for 
flow  to  get  from  the  inflow  boimdary  to  the  outflow  boundary  based  on  the  average  inlet 
velocity.  After  the  jet  had  developed,  the  simulation  was  continued  for  an  additional  5000  time 
steps  (approximately  five  flow-through  times).  KIVA  automatically  adjusts  the  size  of  the  time 
step  by  the  criteria  previously  outlined.  However,  for  the  last  5000  time  steps  it  was  not  allowed 
to  exceed  5.0e-5  seconds. 
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The  instantaneous  flow  flelds  at  two  different  times  are  shown  in  Figures  3.4  and  3.5. 
These  two  flow  fields  are  from  an  axial  plane  of  40x40-cm  cut  from  the  center  of  the  domain  at 
two  different  instances  in  time  0.0 1  seconds  apart.  As  can  be  seen  the  flow  is  highly  unsteady 
and  is  certainly  turbulent.  However,  to  compare  the  numerical  results  obtained  in  this  study  with 
the  experimental  data  requires  that  average  velocity  profiles  be  calculated.  Calculating 
meaningful  averaged  velocity  profiles  for  different  axial  locations  is  not  a  trivial  task.  This  is 
due  to  the  fact  that  the  jet  varies  significantly  in  the  radial  direction.  However,  the  most  difficult 
problem  is  associated  with  the  fact  that  the  Jet  is  non-stationary  in  that  it  fluctuates  about  the 
centerline.  This  is  demonstrated  in  Figure  3  where  the  axial  velocity  at  the  centerline  is  shown  to 
vary  significantly  over  time.  This  is  an  indication  that  the  large-scale  coherent  structures  of  the 
jet  are  highly  unsteady,  and  as  a  result  the  jet  is  rarely  positioned  about  the  center  of  the  domain. 
To  obtain  an  accurate  time  average  many  more  time  steps  would  need  to  be  taken  which  is  not 
feasible  at  this  time.  However,  it  is  believed  that  a  sufficient  number  of  time  steps  have  been 
taken  such  that  a  qualitative  comparison  to  the  experimental  data  can  be  made.  Despite 
oscillatory  nature  of  the  jet,  temporal  and  partially  spatial  averaged  velocity  profiles  were 
calculated  at  two  different  axial  locations  (z  =  20-cm  and  z  =  40-cm).  The  spatial  averaging  was 
performed  by  averaging  the  velocity  profiles  at  four  circumferential  locations  (0  =  0,  0  =  nil, 
0  =  71, 0=37t/2).  The  spatially  averaged  profiles  that  resulted  from  this  procedure  were  then  used 
to  obtain  a  time  average. 

An  important  side  note  to  this  investigation  involves  the  computational  performance  of 
the  ALE  method  and  particularly  the  KIVA  code.  To  perform  the  5000  time  steps  discussed 
above  required  approximately  eight  days  using  a  dedicated  533  MHZ  DEC  ALPHA  NT 
workstation.  It  has  been  demonstrated  that  the  CPU  time  required  to  perform  large  simulations 
using  KIVA  scales  either  linearly,  or  close  to  linear  (Gel  et  al.,  1998).  This  appears  to  be  more 
of  a  feature  of  the  ALE  method  than  of  the  KIVA  code,  and  can  be  considered  a  very  good 
feature.  The  KIVA  code  itself  is  not  particularly  efficient.  The  code  needs  to  be  revised  so  that 
it  can  take  advantage  of  modem  computer  architecture,  which  would  include  paralleli2ation. 
One  of  the  major  drawbacks  of  the  ALE  method  is  the  memory  requirement.  On  minimum 
settings  KIVA  requires  approximately  0.8  kb/node.  However,  the  KIVA  code  was  not  written 
with  LES  in  mind.  It  is  suspected  that  a  code  based  on  the  ALE  methodology  could  be  written 
that  is  more  efficient  with  respect  to  memory  usage. 
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3.4.2  Comparison  to  Experimental  Data 

The  velocity  profiles  obtained  from  averaging  the  numerical  results  are  compared  to  the 
experimental  profiles  in  Figures  3.9-3.20.  From  the  outset  of  this  investigation  it  was  known  that 
the  swirling  jet  problem  was  particularly  sensitive  to  initial  conditions.  This  has  been 
demonstrated  by  Lewellen  et  al.  (1998).  As  a  result  achieving  very  close  agreement  with  the 
experimental  data  should  not  be  expected.  However,  as  can  be  seen  the  overall  the  mean 
velocity  profiles  (Figures  3.9  -  3.14)  compare  fairly  well.  In  particular  the  extent  of  the  radial 
spreading  of  the  jet  is  predicted.  The  simulated  jet  exhibits  a  behavior  that  more  closely 
resembles  circular  jet  farther  downstream  (Figure  3.10).  However,  the  time  average  given  here 
is  somewhat  misleading  due  to  the  oscillatory  nature  of  the  jet  discussed  above.  The  radial 
velocity  is  a  measure  of  the  jet  entrainment.  The  average  radial  velocity  is  given  in  Figures  3.11 
and  3.12  for  z  =  20-cm  and  z=40-cm,  respectively.  It  is  seen  that  the  radial  velocity  is  over¬ 
predicted  at  z=20-cm,  and  this  trend  seems  to  have  been  reversed  at  z=40-cm.  It  is  possible  that 
more  grid  resolution  in  the  radial  plane  may  improve  this  situation.  The  improvement  would  be 
due  to  better  resolution  of  the  shear  region.  More  grid  resolution  in  the  radial  plane  may  also 
improve  the  axial  velocity  agreement  since  the  inflow  velocity  profiles  and  the  resulting  flow 
rate  would  be  better  approximated.  It  is  likely  that  the  simulations  would  then  more  resemble  the 
experimental  situation.  The  swirl  velocity  agrees  fairly  well  at  both  axial  locations  (Figures  3.13 
and  3.14). 

The  averaging  process  previously  described  was  used  to  obtain  statistics  of  resolved 
turbulence.  These  are  compared  to  the  experimental  profiles  shown  in  Figures  3.15-3.20.  In 
general  the  turbulent  statistics  are  imder-predicted.  This  indicates  that  the  subgrid  contributions 
at  this  grid  resolution  are  significantly  large.  However,  given  the  fact  that  it  is  not  possible  to 
initialize  the  inlet  turbulence  such  that  it  is  identical  to  the  experiments,  then  the  agreement  can 
be  considered  quite  good.  Also,  the  initialized  turbulence  was  steady  as  opposed  to  the  real 
situation  where  it  is  unsteady.  To  create  a  more  realistic  simulation,  the  initialized  turbulence 
should  be  both  spatially  and  temporally  correlated.  Another  problem,  which  has  been  previously 
discussed,  is  related  to  the  averaging  process  and  non-stationary  nature  of  the  jet.  The  RMS 
fluctuating  axial  velocity  is  shown  in  Figures  3.15  and  3.16  for  the  two  different  axial  positions. 
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The  levels  are  under-predicted  in  some  regions.  Similar  results  were  obtained  for  the  other  RMS 
fluctuating  velocity  components  (Figures  3.17-3.20). 

3.4.3  Comparison  to  WVU  LES  Results 

As  was  mentioned  previously,  work  on  this  particular  flow  situation  has  been  performed 
using  the  WVU-LES  code.  However,  both  the  grid  and  the  subgrid  scale  model  used  were 
different.  Given  these  differences  exact  duplication  of  the  WVU-LES  results  using  KIVA  is  not 
possible.  However,  qualitative  duplication  of  the  WVU-LES  results  using  KIVA  is  desired  and 
would  yield  further  validation  of  the  KIVA  code. 

Figures  3.6  and  3.7  show  the  instantaneous  pressure  contours  on  the  radial  plane  at  two 
different  axial  locations  (z  =  10-cm  and  z  =  40-cm,  respectively).  Again  the  flow  field  appears 
turbulent  and  is  not  centered  about  the  center  of  the  domain.  Figure  6  shows  that  the  eddies  are 
relatively  symmetric  about  the  center  of  the  domain  at  z  =  10-cm.  Further  downstream  at  z  =  40- 
cm  radial  diffusion  of  the  vorticity  has  occurred  (Figure  3.7).  In  addition  the  structure  of  the 
flow  field  is  tri-modal.  A  similar  structure  is  observed  at  z  =  40  cm  using  the  WVU-LES  code 
(Figure  3.8).  In  addition,  the  order  of  the  variances  between  the  experimental  data  and  the 
simulations  reported  by  Lewellen  et  al.  (1998)  are  quite  similar  to  the  present  results. 

Due  to  the  similarity  between  the  results  obtained  with  KIVA  and  those  obtained  with  the 
WVU-LES  code  it  is  believed  that  some  success  has  been  achieved  in  capturing  the  turbulent 
nature  of  this  flow  field  despite  the  comparatively  low  grid  resolution. 

3.4.4  Validation  of  the  Second-Order  Time  Integration  Scheme 

The  jet  was  simulated  for  an  additional  5000  time  steps  using  the  newly  implemented 
Adams  Bashforth/Runge  Kutta  time  integration  scheme  (see  Chapter  4  for  details).  Averaged 
profiles  were  obtained  in  the  same  manner  as  those  obtained  with  the  first-order  scheme.  Since  it 
is  the  nature  of  the  unsteadiness  that  is  of  interest,  only  the  fluctuating  quantities  obtained  firom 
both  schemes  are  compared.  Also,  since  the  flow  field  is  more  turbulent  closer  to  the  jet  annulus, 
the  closer  of  the  two  axial  locations  will  be  considered  (i.e.  Z=20  cm).  The  three  fluctuating 
components  of  velocity  obtained  using  the  two  different  time  integration  schemes  are  shown  in 
Figures  3.15  to  3.18.  As  can  be  seen  the  profiles  are  nearly  the  same,  with  small  differences 
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being  attributed  to  the  unsteadiness  of  the  flow  field,  not  the  numerics.  This  result  can  be 
justified  by  an  order  of  magnitude  analysis  the  purpose  of  which  is  to  determine  whether  the 
temporal  or  spatial  discretization  errors  are  dominating  the  overall  error.  The  error  assessment 
technique  known  as  Richardson  extrapolation  will  be  used  for  this  purpose.  The  dominating 
error  for  the  temporal  and  spatial  errors  using  the  first  order  backward  Euler  time  integration 
scheme  (i.e.  the  default  in  KIVA)  and  central  differencing  for  the  convective  terms  is  0(At,  Ax^). 

U  =  C,  A/  +  C,  Ax^  «  ^A/  + 

r 


^Ax^ 


(3.2) 


Where  Uref  is  the  reference  velocity,  x  is  the  reference  time  scale,  and  /  is  the  reference 
spatial  scale.  The  reference  time  is  determined  from 

r  =  (3.3) 

which  can  be  substituted  into  the  Equation  3.2  above  to  yield 


U 


r^A/  +  ^Ax' 

««/  I 


The  ratio  of  the  errors  then  becomes 


(3.4) 


temporal  error  ^  ^ 

spatial  error  Ax^ 

The  reference  velocity  is  related  to  the  turbulent  kinetic  energy,  k,  as 

(3.5) 

The  reference  length  scale  for  a  round  jet  is  given  by  Wilcox  (1993)  as 

/  =  0.1  (3.6) 

where  5  is  the  jet  half  width  (jet  radius).  The  time  step  size  and  spatial  discretization  as  well  as 
other  parameters  are  listed  in  Table  3.1  for  different  axial  locations. 

Using  the  experimental  values  for  k  and  5,  Ct/Cx  =  0.13  at  Z  =  20  cm  and  Ct/Cx  =  0.04  at 
Z  =  40  cm.  When  it  is  considered  that  the  sub-cycling  time  step  size  is  1/10  that  of  the  global 
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time  step  given  above  these  nixmbers  drop  by  a  factor  of  ten.  Thus  this  analysis  shows  that  for 
the  chosen  grid  and  time  step  size  the  spatial  error  is  dominant. 


Table  3.1 :  Parameters  needed  for  error  analysis. 


Parameter 

Z  =  20  cm  (2  jet  diameters 

downstream) 

Z  =  40  cm  (4  jet  diameters 

downstream) 

kmax  (m^/  S'^  ) 

experimental 

27.90 

7.52 

numerical 

8.50 

2.57 

6(m)  (experimental) 

0.25 

0.30 

At(s)  (numerical) 

5.0x10'' 

5.0x10'' 

Ax  (m)  (numerical) 

0.007 

0.01 

3.5  Discussions 

This  work  outlines  the  preliminary  attempt  at  LES  of  a  swirling  jet  using  an  ALE 
approach  in  conjunction  with  the  KJVA-3  code.  The  goals  were  to  capture  some  of  the  non¬ 
linearity  and  hence  some  of  the  larger  coherent  turbulent  structures  and  unsteady  effects.  If 
successful  then  some  proof  of  the  LES  capabilities  of  both  the  ALE  method  and  the  KIVA  code 
will  have  been  obtained.  This  is  especially  the  case  given  the  relatively  coarse  grid  used  in  this 
study. 

The  preliminary  results  show  that  the  mean  velocity  profiles  were  fairly  well  predicted. 
Also,  the  simulation  captures  the  unsteady  nature  of  the  jet,  and  predicts  some  of  the  turbulent 
statistics  well.  This  was  achieved  despite  the  fact  that  the  chosen  problem  is  extremely  sensitive 
to  initial  conditions,  and  duplicating  the  experiments  exactly  is  nearly  impossible.  This  is 
certainly  evidence  that  the  both  the  ALE  method  and  the  KIVA  code  in  particular  are  accurate 
enough  to  be  used  for  LES  given  proper  selection  of  convective  differencing  scheme  and  grid. 
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Also,  the  method  does  have  the  potential  for  good  scaling  with  respect  to  grid  resolution.  This 
later  point  makes  the  ALE  approach  attractive  compared  to  other  finite  volume  methods. 

In  addition  to  validating  the  LES  potential  of  the  KIVA  code,  a  new  time  integration 
scheme  has  been  tested.  Time  accuracy  has  been  shown  to  be  relatively  unimportant  for  this 
particular  flow  situation  and  in  particular  with  the  grid  and  time  step  size  used.  However,  time 
accuracy  will  be  important  for  other  problems  and  even  this  problem  when  the  grid  is  refined.  It 
was  demonstrated  that  the  new  second-order  scheme  that  has  been  implemented  into  KIVA  is 
working. 
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Figure  3.1  Grid  in  the  radial  plane. 
Outflow  boundary 


Inflow  boundary 


Figure  3.2  Grid  in  the  axial  plane. 
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;ure  3.9  Average  axial  velocity  at  z=20 


Figure  3.1 1  Average  radial  velocity  at  Z=20 
cm 


Radial  Distance  (cm)  Radial  Distance  (cm) 


Figure  3.10  Average  axial  velocity  at  z=40 


Figure  3.12  Average  radial  velocity  at  z=40 


Figure  3.13  Average  swirl  velocity  at  z=20  Figure  3.15  at  z  =  20  cm. 

cm. 


Radial  Distance  (cm)  Radial  Distance  (cm) 


Figure  3.14  Average  swirl  velocity  at  z=40  Figure  3.16  at  z  =  40  cm. 


Chapter  4  Accuracy  and  Efficiency  Improvements  in  KIVA 


4.1  Accuracy  Improvements 

The  numerical  scheme  in  the  original  KIVA-3  was  kept  simple  and  robust  as  much  as 
possible  in  order  to  incorporate  all  the  capability  to  model  complex  flows  and  make  KIVA-3  a 
general  purpose  code  (i.e.,  laminar/turbulent  flows,  single-phase  or  multiphase  with  fuel  spray 
dynamics,  combustion  chemical  reactions,  etc.).  This  was  done  by  introducing  a  numerical 
scheme  that  is  first  order  accurate  in  time  and  approximately  second  order  accurate  in  space. 
However,  today’s  challenging  problems  demand  more  accurate  schemes  with  higher  resolution 
and  faster  convergence.  Particularly  for  LES  higher  accuracy  of  at  least  second  order  is 
necessary. 

In  this  section  the  improvements  to  enhance  the  spatial  and  temporal  accuracy  of  the 
KIVA  code  are  discussed.  First,  time  integration  scheme  improvements  are  presented.  These 
improvements  were  as  a  result  of  joint  effort  with  Mr.  M.  Prinkey  (1999).  The  second  part  is 
dedicated  to  the  spatial  accuracy  improvements  in  the  code.  A  new  flag  has  been  introduced  to 
run  the  code  with  central  difference  scheme  for  the  convective  terms  in  the  momentum  equations 
while  applying  Quasi-Second  Order  Upwind  (QSOU)  scheme  for  the  density,  energy  and  other 
scalar  equations.  Both  of  these  implementations  are  verified  by  comparing  the  results  with  the 
two-dimensional  lid  driven  cavity  flow  calculations  of  Ghia  et  al.  (1982).  Furthermore  to  assess 
the  overall  spatial  accuracy  of  the  code,  the  simulation  of  a  three-dimensional  lid-driven  cavity 
was  performed  and  compared  with  benchmark  results.  The  details  of  these  results  can  be  found 
inGeletal.(1998). 

4.1.1  Time  Integration  Improvements 

Time  advancement  in  KIVA-3  is  performed  essentially  in  three  phases  (i.e..  Phase  A,  B, 
and  C).  In  Phase  A,  spray  droplet  collision  and  oscillation/breakup  terms  are  calculated  and 
explicit  source  terms  due  to  chemical  reactions,  fuel  spray  etc.  are  integrated  in  the  mass  and 
energy  equations.  Phase  B  calculates,  in  a  coupled,  implicit  fashion,  the  acoustic  mode  terms 
(namely  the  pressure  gradient  in  the  momentum  equation  and  velocity  dilatation  terms  in  mass  & 
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energy  equations),  the  spray  momentum  source  term,  the  source  terms  in  the  turbulence-model 
equations,  and  the  terms  due  to  diffusion  of  mass,  momentum  and  energy.  Phase  C  is  the  rezone 
phase,  in  which  the  convective  transport  associated  with  the  moving  mesh  relative  to  the  fluid  is 
calculated  by  freezing  the  flow  field  and  then  rezoning  (remapping)  it  onto  a  new  computational 
mesh  (in  Lagrangian  phase)  or  onto  the  same  computational  mesh  (in  Eulerian  phase).  All  these 
phases  apply  to  the  same  time  interval  but  one  at  a  time  built  on  top  of  each  other.  The  advective 
terms  of  Phase  C  are  advanced  explicitly  in  consecutive  sub-cycle  timesteps  (Ate)  which  is 
simply  an  integral  submultiple  of  the  global  computational  timestep  (At).  However,  the  diffusion 
terms  of  Phase  B  can  be  integrated  explicitly,  implicitly,  or  semi-implicitly.  In  certain  problems, 
when  timesteps  are  small  enough,  KIVA-3  automatically  uses  a  stable  explicit  scheme  for  which 
no  costly  iterative  solution  is  required.  When  timesteps  are  not  suitable  for  explicit  marching,  a 
semi-implicit  time  advancement  option  is  employed.  The  degree  to  which  the  diffusion  terms 
are  implicitly  discretized  is  based  on  a  combination  of  stability  and  efficiency  considerations. 
The  implicitness  parameter  (()>)  enables  some  sort  of  weighting  of  the  old-  and  new-time  values 
of  the  solution  variable. 

Numerical  experiments  have  shown  that  this  approach  improves  the  computational 
efficiency  without  degrading  the  numerical  stability  when  possible.  Two  separate  implicitness 
factors  are  calculated  in  KIVA-3.  The  first  one  is  (j)d  which  is  associated  with  the  diffusion 
terms.  The  second  one  is  <j)p  which  is  associated  with  the  pressure  wave  propagation. 

During  Phase  C,  sub-cycle  timesteps  are  taken  to  advance  the  convective  terms  explicitly. 
The  convective  terms  are  advanced  with  time  steps  that  are  the  same  size  or  smaller  than  what  is 
used  for  the  diffusion  terms.  The  time  step  size  of  the  each  sub-cycle  is  based  on  stability 
considerations.  Though  the  overall  time  accuracy  is  only  first-order  accurate,  the  global  time 
step  is  based  on  several  considerations  including  stability  and  accuracy  constraints  which 
restricts  the  time  step  to  fairly  small  values  hence  assuring  accuracy. 

Further  details  on  the  time  integration  scheme  and  the  implicitness  parameters  can  be 
found  in  Appendix  H  of  the  KIVA-II  manual  (Amsden  et  al.  (1989)). 
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The  new  second  order  accurate  time  scheme 


For  the  problems  without  any  spray/droplet  and  chemical  reactions,  the  source  term 
contributions  due  to  these  terms  can  be  omitted.  Also  if  we  consider  that  there  are  no  source 
terms  in  Phase  B  this  yields  the  semi-implicit  formulation  for  the  diffusion  terms.  Following  the 
suggestion  in  Appendix  H  of  KIVA-II  manual  (Amsden  et  al.  (1989)),  setting  the  implicitness 
factor  for  diffusion,  (|)d ,  equal  to  0.5  throughout  the  simulation  yields  to  second  order  accurate 
Crank-Nicolson  scheme  which  may  have  overshoots  and  lead  to  inaccurate  results  for  large 
values  of  Cd  (=  vAt  /  Ax^).  Hence  it  follows  that  by  selecting  an  appropriate  valued  for  <t)d,  the 
KIVA-3  code  can  be  made  second  order  accurate  in  time  except  only  the  time  integration  steps  in 
Phase  C. 

For  this  purpose  a  scheme  consisting  of  a  two  stage  Runge-Kutta  and  Adams-Bashforth 
method  was  proposed  and  implemented  for  the  advected  quantities  (e.g  internal  energy,  pE  and 
momentum,  pU)  in  Phase  C.  Further  details  of  the  implementation  of  the  new  time  scheme  can 
be  found  in  Gel  (1999). 

4.1 .2  Spatial  Accuracy  Improvements 

All  derivatives  other  than  the  convective  terms  are  approximated  with  central  differences. 
Spatial  accuracy  of  the  convective  terms  is  limited  to  several  choices  including  first-order 
upwinding,  central  differencing,  a  user  defined  scheme  that  combines  the  previous  two  choices 
(PDC),  and  a  flux-limited  higher-order  scheme  known  as  Quasi-Second-Order  Upwinding 
(QSOU).  Both  the  first-order  upwinding  and  the  QSOU  schemes  are  monotonic.  It  is  important 
to  note  that  KIVA-3  is  a  compressible  flow  code.  Oscillations  from  non-monotonic  schemes  can 
be  catastrophic  for  compressible  codes  where  scalar  variables  such  as  density  and  temperature 
are  involved.  This  can  lead  to  code  instability  and/or  largely  inaccurate  solutions.  As  a  result, 
monotonic  schemes  are  the  preferred  choice  in  KIVA-3.  The  trade-off  for  monotonacity  is  the 
loss  of  accuracy  near  steep  gradients.  This  is  true  for  all  monotonic  schemes,  including  what  is 
probably  the  most  accurate  of  these  type  of  schemes  known  as  the  Piecewise  Parabolic  Method 
(PPM)  (Colella  and  Woodward,  1984).  On  a  uniform  grid  in  the  absence  of  steep  gradients,  PPM 
is  formally  fourth-order  accurate  whereas  the  QSOU  method  is  second-order  accurate.  In  the 
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presence  of  steep  gradients  PPM  reduces  to  second-order  accuracy  whereas  QSOU  reduces 
nominally  to  a  kind  of  first-order  upwind  scheme. 

Another  convective  scheme  option  included  in  the  Partial  Donor  Cell  (PDC)  formulation 
in  KIVA-3  is  the  central  differencing  (CD).  Applying  CD  is  difficult  due  to  the  well-known 
dispersive  error  that  results  when  this  scheme  is  applied  to  a  situation  where  the  cell  Peclet 
number  is  greater  than  2.0.  This  constraint  imposes  prohibitively  fine  grids.  When  CD  is  used 
special  care  needs  to  be  taken  to  prevent  numerical  dispersion  from  dominating  the  solution. 
This  includes  using  some  form  of  deferred  correction,  or  explicitly  letting  extra  diffusion  enter 
into  the  calculation.  In  practice  higher-order  upwinding  methods  (such  as  the  QUICK  scheme) 
are  often  used  to  obtain  grid  independent  results  with  a  much  smaller  grid.  However,  no  such 
scheme  is  available  in  KIVA-3,  and  for  the  reasons  that  will  be  discussed  in  the  next  paragraph 
upwinding  schemes  are  not  specifically  desirable. 

CD  is  of  particular  interest  since  it  is  typically  used  when  finite  volume  methods  are 
applied  in  the  Large  Eddy  Simulation  (LES)  approach  (Piomelli,  1998;  Mittal  and  Moin,  1997). 
CD  is  used  in  order  to  prevent  numerical  diffusion  from  polluting  the  calculations.  Truncation 
errors  associated  with  numerical  diffusion  are  more  harmful  to  LES  than  those  that  are  associated 
with  numerical  dispersion.  Unfortunately,  the  dominating  errors  that  result  from  any  kind  of 
upwinded  scheme  are  diffusive.  The  fine  grid  resolution  used  in  LES  combined  with  the  added 
diffusion  provided  by  the  subgrid-scale  turbulence  model  causes  the  effective  Peclet  number  to 
be  smaller  than  it  would  be  in  the  laminar  case.  This  allows  some  control  over  the  dispersion 
errors  without  the  introduction  of  subgrid-scale  damping  numerical  viscosity. 

KIVA-3  has  a  fourth-order  node  coupler  that  is  included  as  part  of  the  base  code.  This 
routine  is  provided  as  a  means  of  filtering  oscillations  that  occur  due  to  grid-to-grid  de-coupling 
of  the  pressure  and  the  velocity  fields.  Pressure- velocity  de-coupling  occurs  when  CD  is  used  in 
the  pressure  solver  on  a  non-staggered  grid.  The  current  release  of  KIVA-3  should  not  encounter 
the  pressure  de-coupling  problems  that  existed  with  the  older  versions  since  a  staggered  grid  is 
now  used.  However,  the  node-coupler  is  still  included  in  the  code.  This  routine  appears  to 
introduce  fourth-order  diffusion,  which  is  used  to  filter  the  oscillations.  It  should  be  recognized 
that  the  oscillations  that  result  from  pressure  de-coupling  are  similar  to  those  that  result  from 
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applying  CD  to  the  convective  terms.  The  use  of  the  alternate  node  coupler  as  a  filter  to  remove 
those  unphysical  oscillations  that  result  from  using  CD  is  proposed  as  an  alternative. 

Several  parametric  studies  to  assess  the  degree  to  which  the  node-coupler  needs  to  be 
employed  to  remove  the  unphysical  oscillations  were  carried  out.  In  order  to  assess  the  accuracy 
of  the  overall  scheme,  a  lid-driven  cavity  test  case  was  simulated  with  KIVA-3  after  some 
modifications  in  the  code. 

For  an  LES  application  the  added  turbulent  eddy  viscosity  is  usually  enough  to  ensure 
that  the  effective  Peclet  numbers  are  small  enough  to  minimize  the  dispersive  errors  in  the 
momentum  equations  resulting  from  CD. 

Small  oscillations  in  density  and  energy  can  grow  and  lead  to  non-convergence  of  the 
solution.  However,  for  truly  incompressible  problems  such  as  that  studied  here  (i.e.,  lid  driven 
cavity  problem)  the  energy  equation  is  de-coupled  firom  the  velocity  field  and  both  density  and 
temperature  should  remain  constant.  To  ensure  that  this  was  the  case  a  modification  was  made 
to  the  code.  Whenever  CD  is  used  for  the  convective  terms  in  the  momentum  equation,  QSOU  is 
used  for  all  scalar  variables.  This  appears  to  have  negligible  effect  on  the  simulations  since  both 
temperature  and  density  remain  constant,  except  for  very  small  pertxirbations  (4  %  difference 
between  maximum  and  minimum  density  for  the  problem  studied  here).  This  is  as  would  be 
expected  since  these  quantities  should  be  de-coupled  from  the  velocity  field  and  should  not 
change. 

4.1.3  Verification 

The  time  and  spatial  accuracy  improvements  presented  mentioned  above  were  verified  in 
case  of  the  well  known  benchmark  problem  of  a  lid  driven  cavity.  The  results  of  this  verification 
work  can  be  found  in  Gel  (1999).  A  separate  study  to  assess  the  overall  spatial  accuracy  of  the 
code  with  the  spatial  accuracy  improvement  was  performed  with  the  three-dimensional  lid  driven 
cavity  test  case.  The  results  can  be  found  in  Gel  et  al.  (1998). 

In  addition  to  the  above  test  cases,  the  time  and  spatial  accuracy  improvements  presented 
in  this  study  were  utilized  and  used  in  the  LES  of  a  free  swirling  jet  in  the  study  of  Smith  et  al. 
(1999a).  Although  numerical  accuracy  improvement  analysis  was  not  conducted,  the  fact  that  the 
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new  implementations  gave  very  close  results  or  sometimes  better  results  than  the  original, 
indicates  the  validity  of  the  modifications  introduced. 

4.2  Efficiency  Improvements 

In  an  attempt  to  improve  the  computational  efficiency  of  KIVA-3,  the  bottlenecks  in  the 
single  processor  (i.e.,  original)  version  were  investigated  before  proceeding  with  the  parallel 
version  implementation.  For  this  purpose  the  two  example  cases  provided  with  the  KIVA-3 
manual  were  evaluated,  performance  wise,  at  the  subroutine  level  profiling.  These  particular  test 
cases  were  selected  due  to  the  fact  that  they  constitute  the  majority  of  the  features  of  the  code 
during  the  execution.  As  seen  from  the  Table  4.1,  implicit  pressure  solver  (PSOLVE,  RESP)  is 
one  of  the  most  CPU  intensive  segments  of  the  code  which  becomes  more  crucial  in  pressure 
dependent  problems.  The  last  column  in  this  table  shows  about  13  %  of  the  total  execution  time 
is  spent  on  that  particular  subroutine. 

Table  4. 1  Listing  of  top  ten  routines  with  the  highest  CPU  utilization  during  the  execution  of 
TEAPOT  example  case 


Subroutine 

Tot  Time 

Calls 

Avg.  Time 

% 

CCFLUX 

2.58e+01 

1649 

1.56e-02 

28.75 

RESP 

1.09e+01 

4257 

2.55e-03 

12.11 

MOMFLX 

8.69e+00 

1649 

5.27e-03 

9.68 

REST 

4.85e+00 

1778 

2.73e-03 

5.40 

PSOLVE 

4.48e+00 

119 

3.76e-02 

4.99 

TSOLVE 

3.35e+00 

119 

2.81e-02 

3.73 

RESUVW 

2.67e+00 

532 

5.03e-03 

2.98 

EXDI 

2.48e+00 

no 

2.25e-02 

2.76 

BC 

2.14e+00 

2981 

7.19e-04 

2.39 

RESK 

2.09e+00 

1311 

1.59e-03 

2.32 

An  extensive  investigation  of  the  implicit  iterative  solvers  in  KIVA-3  was  carried  out  due 
to  lack  of  concise  documentation  on  the  implementation  of  Conjugate  Residual  (CR)  method  in 
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the  code.  Improving  the  preconditioner  used  in  the  original  version  was  determined  to  have  a 
significant  impact  on  single  processor  performance. 

In  this  section  the  implementation  of  a  more  sophisticated  preconditioner  and  the 
consequent  performance  improvement  without  any  overhead  in  memory  allocation  is  discussed. 
The  implementation  was  carried  out  in  such  a  way  that  the  user  can  simply  replace  the  two 
associated  subroutines  with  the  new  ones  to  get  the  new  preconditioner  implemented  into  their 
version  of  KIVA-3.  The  details  of  the  improved  pressure  solver  can  be  found  in  Smith  et  al. 
(1999b). 

4.2.1  Conjugate  Residual  Method 

KIVA-3  features  a  Conjugate  Residual  (CR)  method  (Amsden  et  al.  (1985))  based 
algorithm  for  the  solution  of  matrices  iteratively.  This  solution  approach  is  activated  only  when 
the  value  of  implicitness  factor  is  greater  than  zero,  i.e.  semi-  or  fully-implicit.  For  example  for 
pressure  ,  if  (jip  =  0.0  then  the  solution  of  the  pressure  field  is  obtained  explicitly.  However,  if  (|»p 
>  0.0  then  the  solution  is  obtained  iteratively  with  the  CR  method  based  solver  (PSOLVE).  This 
method  falls  in  the  class  of  matrix  solution  procedures  known  as  “Krylov  subspace”  methods 
(see  Saad  (1984)  and  Golub  et  al.  (1996)  for  further  information  on  Krylov  subspace  methods). 
Since  the  initial  implementation  of  the  CR  method  in  KIVA  in  1986,  significantly  improved 
methods  like  Conjugate  Gradient  and  Multigrid  methods  have  been  developed  and  widely  used. 
However,  CR  method  was  not  replaced  with  any  of  these  newer  methods  due  to  the  memory 
conservative  approach  followed  in  implementing  it  into  KIVA,  i.e.,  the  coefficient  matrix  is  not 
stored. 

Another  reason  that  has  hindered  the  implementation  of  improved  algorithms  is  probably 
the  lack  of  detailed  documentation  on  the  solver  algorithms  of  KIVA-3.  This  shortcoming  has 
necessitated  considerable  effort  to  understand  the  implementation  logic  of  the  CR  method  in 
order  to  be  able  to  modify  the  code  to  improve  its  efficiency.  A  detailed  description  of  the 
pressure  solution  technique  employed  in  KIVA-3  is  given  in  Gel  (1999)  and  Smith  et  al.  (1999b). 
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4.2.2  Symmetric  Gauss-Seidel  and  Successive  Over-Relaxation  Preconditioning 
(SGS-SSOR) 

The  traditional  way  of  generating  a  precondition  matrix  is  to  find  an  approximate  inverse 
of  the  coefficient  matrix.  This  could  be  done  in  several  ways.  The  simplest  preconditioning 
method  is  known  as  “Jacobi  preconditioning”,  which  is  also  used  in  KIVA-3  but  for  reasons  that 
will  become  apparent  later,  this  method  is  more  appropriately  referred  to  as  “diagonal 
preconditioning”.  In  this  method,  the  reciprocal  of  the  center  diagonal  of  coefficient  matrix.  A, 
is  used  as  preconditioner.  This  approach  is  memory  efficient  since  only  one  full  size  array  needs 
to  be  saved.  However,  it  is  a  primitive  form  of  preconditioning.  More  advanced  procedures 
involve  calculating  an  incomplete  LU  (ILU)  decomposition,  or  an  incomplete  Cholesky  (IC) 
factorization  of  the  original  matrix.  The  later  is  a  special  case  of  ILU  decomposition  for 
symmetric  matrices.  The  coefficient  matrix  is  approximately  factored  such  that  only  non-zero 
elements  that  exist  in  the  original  coefficient  matrix  are  considered.  Although  this  procedure 
often  yields  very  good  preconditioning,  the  generation  of  the  incomplete  factorization  is  CPU 
intensive  and  is  not  easily  parallelizable.  Furthermore,  the  memory  allocation  is  increased  due  to 
fact  that  the  factored  matrix  needs  to  be  stored.  For  the  pressure  matrix  in  KIVA-3  this  would 
require  an  additional  seven  full  size  arrays  which  is  a  significant  burden  for  a  memory  bounded 
code.  A  substantially  improved  approach  for  preconditioning  can  be  implemented  by  performing 
several  sweeps  of  a  simple  iterative  solution  procedure  such  as  Jacobi,  Gauss  Siedel  (GS),  or 
Successive  Over-Relaxation  (SOR)  and  in  this  approach  only  the  original  coefficient  matrix  need 
to  be  stored. 

The  advantage  of  the  GS/SOR  approach  is  that  convergence  is  accelerated  over  the 
Jacobi  method.  Also,  only  one  array  is  required  for  the  later  two  methods  since  new  values  are 
simply  substituted  for  the  old  values.  However,  two  arrays  are  required  for  the  Jacobi  approach 
since  two  iteration  levels  must  be  stored. 

In  most  cases  the  pressure  matrix  will  be  non-symmetric  due  to  various  reasons  such  as 
the  non-uniformity  of  grids  used.  However,  for  some  special  cases  the  matrix  will  be  symmetric, 
for  which  the  Krylov  methods  converge  faster.  In  fact,  for  many  Krylov  subspace  methods, 
including  the  CR  method,  convergence  is  guaranteed  only  if  the  matrix  is  symmetric.  This 
seldom  causes  problems  in  application  in  that  the  procedure  rarely  fails  to  find  a  solution. 
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Depending  on  how  it  is  applied,  preconditioning  may  destroy  the  symmetry  of  the  matrix.  One 
directional  sweeps  of  the  Jacobi  method  will  retain  symmetry.  However,  one  directional  sweeps 
of  both  GS  and  SOR  will  destroy  symmetry.  For  those  cases  where  the  matrix  is  symmetric  it  is 
desirable  to  apply  the  preconditioning  in  such  a  way  as  to  preserve  symmetry.  This  can  be  done 
by  using  symmetric  GS  (SGS)  and  symmetric  SOR  (SSOR).  Employing  symmetric  versions  of 
these  schemes  is  a  simple  procedure  that  entails  making  two  separate  sweeps;  the  first  from  the 
first  node  in  the  matrix  to  the  last,  then  reversing  the  sweep  direction  from  the  last  node  to  the 
first.  The  two  opposite  direction  sweeps  represent  one  iteration  of  the  symmetric  scheme.  A 
detailed  discussion  on  the  implementation  issues  of  this  new  combined  approach  and  how  the 
original  memory  allocation  requirements  of  KIVA-3  were  preserved  without  introducing  any 
memory  overhead  can  be  found  in  Gel  (1999). 

4.2.3  Benchmark  Problems 

Four  different  benchmark  problems  were  used  to  demonstrate  the  effectiveness  of  the 
new  preconditioning  scheme.  The  types  of  problems  vary  from  axisymmetric  (two-dimensional) 
to  full  three-dimensional  simulations.  The  test  cases  also  include  compressible  flows  with 
combustion,  which  are  more  relevant  to  engine  applications,  and  a  flow  situation  which  is  mostly 
incompressible. 

The  first  problem,  hereafter  referred  to  as  Case  1,  was  the  standard  Baseline  case  that  is 
included  as  an  example  problem  with  KIVA  distribution.  This  problem  is  the  case  of  a  two- 
dimensional  (axisymmetric)  engine  cylinder  with  compression  and  expansion  strokes  and 
including  combustion.  There  were  a  total  of  1,380  grid  cells  in  this  problem. 

The  second  problem  (Case  2)  is  the  same  as  Case  1,  except  that  the  axisymmetric  grid  has 
been  rotated  360°  to  make  it  three-dimensional.  The  rotation  was  achieved  by  placing  36  grid 
points  in  the  tangential  direction,  giving  approximately  16,000.  As  with  the  first  problem,  only 
the  compression  and  expansion  strokes  are  simulated,  but  combustion  is  still  included.  This 
problem  should  demonstrate  how  the  pressure  solution  dominates  the  overall  solution  time  as 
both  the  number  of  grid  cells  and  the  dimensionality  of  the  problem  increases. 

The  third  problem  (Case  3)  is  the  same  as  Case  2,  but  72  grid  points  are  placed  in  the 
tangential  direction.  Thus  the  total  number  of  grid  points  is  twice  that  of  Case  2  at  approximately 
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32,000.  This  problem  is  used  to  test  the  performance  of  the  new  preconditioning  scheme  as  the 
grid  resolution  becomes  finer.  It  is  expected  that  the  pressure  solver  will  dominate  the  overall 
CPU  time  more  than  it  does  in  Case  2.  As  a  result,  the  quality  of  preconditioning  will  become 
more  important.  This  phenomenon  is  related  to  the  fact  that  the  matrix  becomes  more  difficult  to 
solve  (invert)  as  it  becomes  larger. 

/ 

The  fourth  problem  (Case  4)  is  a  different  kind  of  problem  from  the  first  three  in  that  it 
does  not  involve  an  engine  cylinder.  It  is  a  three-dimensional  swirling  jet  defined  on  a  non- 
uniform  rectangular  grid.  A  total  of  64,000  grid  cells  were  used.  This  particular  problem  does 
not  include  combustion  and  is  for  the  most  part  incompressible.  It  is  suspected  that  the 
computational  effort  will  be  dominated  by  the  pressure  solution.  The  reader  is  referred  to  Smith 
et  al.  (1998)  for  more  details  on  this  particular  problem. 

Table  4.2  Performance  comparison  of  KIVA-3  and  KIVA-3/WVU  for  the  benchmark  problems 


KIVA-3 

K1VA-3/WVU 

Case 

Number 

Normalized 
CPU  Time 

Pressure 

Iterations 

Normalized 

CPU  Time 

Pressure 

Iterations 

Speedup 

1 

1.0 

38 

0.93 

15 

6.8 

2 

1.0 

54 

2  sweeps  0.81 

29 

18.6 

3  sweeps  0.80 

24 

20.4 

3 

1.0 

75 

2  sweeps  0.90 

51 

9.8 

3  sweeps  0.84 

44 

16.1 

4  sweeps  0.82 

39 

18.3 

4 

1.0 

45 

2  sweeps  0.80 

15 

20.0 

3  sweeps  0.79 

12 

21.3 

The  size  of  the  time  step  in  all  four  problems  is  controlled  by  the  stability /accuracy 
constraints  contained  in  the  original  code.  The  k-e  model  was  used  for  Cases  1-3.  A 
Smagorinsky  turbulence  model  was  employed  in  Case  4,  thus  no  turbulence  transport  equations 
need  to  be  solved  for  this  problem.  All  computations  were  performed  on  a  DEC  Alpha  533  MHz 
workstation  running  Windows  NT  v4.0. 
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4.2.4  Discussion 


The  subroutines  DRDP  and  PSOLVE  have  been  modified  in  such  a  way  that  the  original 
routines  in  both  KIVA-3  and  KIVA-3V  can  be  directly  replaced  to  yield  a  more  efficient  code  on 
a  single  processor.  The  improvement  in  the  efficiency  of  the  overall  code  is  problem  dependent. 
For  problems  that  do  not  contain  a  large  number  of  grid  cells  the  diagonally  preconditioned  CR 
routine  is  robust.  However,  as  the  grid  is  refined  and  the  total  number  of  grid  cells  increases,  the 
solution  of  the  matrix  becomes  more  difficult.  Hence,  the  quality  of  the  preconditioning 
becomes  more  important,  and  the  modifications  outlined  in  this  paper  may  provide  a  significant 
reduction  in  overall  CPU  time.  The  PSOLVE  routine  has  been  re-written  in  such  a  way  that  the 
user  can  choose  between  the  original  preconditioning  scheme  and  the  one  presented  in  this 
chapter. 

Table  4.2  shows  the  speedup  achieved  after  the  implementation  of  the  new  preconditioner 
based  on  the  number  of  sweeps  employed.  The  dimensionality  of  the  problem  seems  to  be  the 
most  important  factor  in  the  effectiveness  of  the  new  preconditioning  scheme.  The  overall  code 
speedup  is  significantly  less  for  the  axisymmetric  problem  (Case  1)  compared  to  the  three- 
dimensional  problems  (Cases  2-4).  This  is  likely  due  to  the  fact  that  the  matrix  produced  is 
relatively  small  and  the  existing  diagonal  preconditioning  is  robust.  Also,  the  matrix  produced 
by  an  axisymmetric  problem  is  easier  to  solve  since  there  are  two  less  diagonals.  The  modified 
pressure  solver  appears  to  have  a  significant  impact  on  three-dimensional  problems  where  the 
number  of  grid  cells  is  large  and  the  quality  of  the  preconditioning  becomes  more  important.  In 
this  work  KIVA-3/WVU  was  observed  to  be  up  to  20  %  faster  than  the  original  code  on  two 
different  kinds  of  three-dimensional  problems.  It  should  be  noted  that  the  CR  method  is  also  used 
in  several  other  routines  including  KESOLVE  (for  k  and  e),  TSOLVE  (for  T),  VSOLVE  (for 
velocities),  and  YSOLVE  (for  species  densities  and  mass  fractions).  Implementing  the  new 
preconditioning  method  in  these  routines  is  expected  to  be  a  simpler  procedure.  It  was  observed 
that  TSOLVE  and  BCESOLVE  are  expensive  for  three-dimensional  problems  such  as  Cases  2  and 
3.  The  modifications  outlined  in  this  work  could  provide  significant  CPU  savings  for  problems 
where  the  local  iterations  for  the  scalar  field  variables  dominate  the  overall  solution  procedure. 
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4.3  Performance  Improvements 


The  KIVA  family  of  codes  have  been  primarily  developed  on  Cray  vector 
supercomputers  with  explicit  directives  to  take  advantage  of  vector  processing  units  without 
stalling  the  pipelines.  However,  the  gain  from  vector  supercomputers  are  becoming  obsolete  due 
to  the  significant  improvements  in  processor  chip  design  and  introduction  of  efficient  parallel 
processing  platforms.  As  a  result  of  these  developments,  the  definition  of  “high  performance 
computing”  has  evolved  and  changed  significantly.  These  changes  had  significant  impact  on  the 
problem  complexity  and  size  to  be  tackled. 

Since  the  release  of  it's  first  version  in  1985,  the  KIVA  code  has  been  also  evolving  and 
improving  in  parallel  to  all  the  new  developments  in  computational  resources  and  CFD 
techniques.  For  example,  the  block-structured  grid  capability  in  KIVA-3  provides  significant 
advantages  over  KIVA-II,  especially  for  complex  internal  combustion  engine  geometries. 
Another  example  is  the  indirect  addressing  introduced  in  KIVA-3,  which  eliminates  the 
obligation  to  go  through  the  storage  array  elements  in  a  particular  sequence.  In  spite  of  these 
improvements  in  the  metholodogy,  there  are  certain  bottlenecks,  which  are  quite  difficult  to 
overcome.  For  example  the  ghost  cells  surrounding  the  grid  still  impose  severe  restrictions  on  the 
maximum  grid  resolution  that  can  be  achieved  on  single  processor  systems. 

Considering  today's  challenging  problems  which  require  adequately  fine  grid  resolution 
capability,  the  most  economical  way  to  achieve  progress  towards  this  goal  with  KIVA-3  is  to 
acquire  the  capability  to  run  this  code  on  parallel  systems,  preferably  on  distributed-memory 
architectures. 

A  shared-memory  parallel  implementation  of  KIVA-3,  which  is  particularly  optimized 
for  Silicon  Graphics  symmetric  multiprocessor  systems,  has  been  available  to  public  for  some 
time.  However,  the  speedup  characteristic  of  this  version  is  not  favorable.  Figure  4. 1  shows  the 
speedup  plot  up  to  eight  processors  for  a  test  case  based  on  a  180°  wedge  with  a  grid  of  80  x  18  x 
80.  It  is  seen  that  speed-up  efficiency  stalls  after  four  processors.  Due  to  this  poor  scalability  and 
parallel  efficiency,  a  domain  decomposition  based  approach  with  explicit  message-passing  was 
investigated  by  Yasar  et  al.  (1992). 
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Figure  4.1  Speedup  plot  of  SGI  optimized  shared-memory  implementation  of  KIVA-3 


Figure  4.2  One-dimensional  domain  decomposition  applied  to  a  sector  domain 
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Figure  4.3  Illustration  of  domain  decomposition  over  3  processors  for  single  layer  in  z-direction 


In  this  study,  a  domain  decomposition  strategy  in  the  x-direction  (or  radial  direction)  is 
based  on  the  idea  of  overlapping  the  left  &  right  faces  of  adjacent  subdomains.  This  approach 
was  initially  proposed  by  Yasar  (1998)  who  developed  a  preliminary  implementation  of  KIVA-3 
on  Intel  Paragon  platform  with  NX  message-passing  libraries. 

The  parallel  version  has  been  ported  to  various  distributed-memory  and  shared-memory 
parallel  computer  platforms  with  MPI  message-passing  library.  Note  that  although  the 
fundamental  ideas  of  the  current  implementation  are  based  on  the  previous  work  by  Yasar  et  al. 
(1992),  the  current  work  essentially  originated  from  scratch  by  building  on  top  of  the  original 
source  code  of  KIVA-3.  The  current  domain  decomposition  approach  is  applicable  to  various 
geometries  in  engine  applications  ranging  from  full  360°  mesh  or  a  sector  mesh  as  shown  in 
Figure  4.2. 

Figure  4.3  shows  how  a  cartesian  computational  domain  of  9  cells  in  length,  2  cells  in 
width  and  height  is  decomposed  into  3  subdomains;  here  in  Figure  4.3,  only  a  single  z-layer  is 
shown  for  brevity.  Different  symbols  have  been  used  to  distinguish  various  cell  types  and 
boundary  vertices  available  in  original  KIVA-3.  The  symbols  can  be  interpreted  as  follows, 

•  Square  symbol  implies  GHOST  cells  on  the  left  and  bottom  faces  of  the  whole 
domain  with  F  =  0  and  FV  =  0  attributes. 

•  Circle  symbol  implies  ACTIVE  cells  with  F  >  0  and  FV  >  0  attributes. 

•  Triangle  symbol  implies  the  BOUNDARY  vertices  (i.e.,  active  vertices  without 
any  cell  association)  with  F  =  0  and  FV  >  0  attributes. 

Domain  decomposition  is  carried  along  the  x-direction  and  each  subdomain  contains  3  x 
2x2  cells.  Each  subdomain  is  assigned  to  a  processor.  Left  and  right  faces  of  each  subdomain 
overlap  with  the  adjacent  subdomains  right  and  left  face.  Each  processor  has  it's  own  grid  vvith 
defined  boundaries  and  runs  the  same  KIVA-3  code  separately.  Note  that,  in  practice  when  a 
computational  domain  is  decomposed  over  several  processors,  grid  resolution  is  increased  as 
high  as  possible  to  take  advantage  of  each  processor.  In  this  example  overall  grid  resolution  was 
kept  low  for  demonstration  purposes  only.  An  enlarged  view  of  the  overlapping  face  between 
adjacent  subdomains  (i.e.,  on  Processors  0  and  1)  is  shown  in  Figure  4.4.  The  vertical  arrows 
between  subdomains  illustrate  the  left/right  faces  where  two  way  communication  between 
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processors  take  place.  KIVA-3  employs  staggered  grids,  i.e.,  scalar  quantities  are  at  cell  centers 
and  vectorial  quantities  (U,  V,  W)  are  at  the  vertices.  This  configuration  obliged  the  necessity  of 
two  distinct  communication  types  to  be  created.  The  first  one  was  for  the  cell-center  based 
quantities  and  the  second  for  the  vertex  based  quantities.  A  more  detailed  illustration  of  various 
Communication  types  among  different  processors  can  be  found  in  Gel  (1999). 

After  completing  the  insertion  of  all  communication  directives  into  the  code,  a  rigorous 
debugging  of  the  parallel  version  was  necessary  to  make  sure  that  all  of  the  necessary  variables 
are  communicated  at  the  appropriate  locations  and  all  the  processors  were  behaving  in  a 
synchronized  manner  without  any  problem. 

4.3.1  Benchmarking  for  Speedup 

An  important  component  of  effective  parallel  computing  is  determining  whether  the 
program  is  running  efficiently  (i.e.,  how  effectively  is  the  processor  being  utilized),  and  how  well 
it  will  perform  if  a  large  number  of  processors  are  used  or  how  well  it  scales  up.  These  issues  are 
usually  quantified  by  two  important  parameters  in  parallel  computing. 

Speedup  =  SerTime(n)  /  ParTime(n,p) 

Efficiency  =  SerTime(n)  /  (p  •  ParTime(n,p)) 

where  SerTime(n)  is  the  time  of  best  serial  algorithm  to  solve  instance  of  a  problem  A  of  the 
size  n  on  one  processor.  ParTime(n,p)  is  the  time  needed  by  a  given  parallel  algorithm  and  given 
parallel  architecture  to  solve  an  instance  of  problem  A  of  size  n,  using  p  processors. 

Note  that  SerTime(n)  <  ParTime(n,l)  and  in  general,  one  expects  to  have 

0<  Speedup  <  p 

0  <  Efficiency  <1.0 

In  order  to  determine  the  efficiency  of  the  distributed-memory  implementation  of  KIVA- 
3  developed  in  this  study,  a  series  of  benchmark  runs  with  different  size  grids  and  total  number 
of  processors  were  carried  out.  All  the  runs  presented  in  this  section  were  conducted  on  a  112 
processor  SGI  Origin  2000  system  at  Department  of  Defense's  ERDC  (formerly  CEWES)  Major 
Shared  Resource  Center  at  Vicksburg,  MS.  The  runs  performed  for  benchmarking  purposes  were 
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conducted  for  a  pre-determined  number  of  timesteps  where  the  dependency  of  speedup  trends  on 
the  number  of  timesteps  were  also  investigated.  It  was  determined  that  beyond  100  timesteps  no 
significant  difference  in  speedup  trends  was  observed.  Due  to  the  significant  overhead  in  setup 
phase  in  single  processor  runs,  speedup  and  parallel  efficiency  figures  were  calculated  first  by 
including  and  then  excluding  the  time  elapsed  in  setup  phase.  Such  a  distinction  was  necessary 
due  to  the  superlinear  speedup  behavior  observed  when  setup  phase  time  was  included. 

The  first  benchmark  problem  is  the  three-dimensional  laminar  channel  flow  at  a  Reynolds 
number  Re=  1,000  which  was  tested  up  to  48  processors  with  several  grid  resolutions  starting 
from  250,000  vertices  to  1,500,000  vertices,  referred  as  250K  and  1500K  grid,  respectively.  This 
problem  was  particularly  selected  for  debugging  purposes  during  the  development  of  the  parallel 
version  and  utilizes  nearly  80  %  of  the  subroutines,  which  was  expected  to  give  an  idea  about  the 
overall  performance.  The  results  show  that  under  favorable  conditions  with  a  simple  problem 
like  laminar  channel  flow,  a  speedup  factor  of  38  could  be  achieved  with  48  processors  for  the 
1500K  grid  case  (see  Figure  4.5)  which  corresponds  to  a  parallel  efficiency  of  0.79.  Note  that  the 
ratio  of  computational  volume  to  communicated  subdomain  face  area  is  an  important  factor  in 
determining  the  optimal  number  of  processors  for  a  given  grid  size. 


Speedup  Plot  for  1500K  Grid 
(excluding  setup  phase  time) 


Figure  4.5  Speedup  plot  for  1500K  cases  based  on  excluding  setup  phase  time. 
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The  second  problem  was  the  turbulent  round  jet  problem  with  an  inlet  velocity  of  1500 
cm/s  and  a  diameter  of  10  cm  enters  a  three-dimensional  rectangular  domain  through  the  left 
wall.  This  problem  is  the  simplified  version  of  the  free  swirling  annular  jet  case  (i.e.,  annular  jet 
replaced  with  round  jet  without  any  swirl)  studied  by  Smith  et  al.  (1999a)  on  a  single  processor 
DEC  Alpha  machine  with  KIVA-3  up  to  500,000  vertices.  Benchmark  tests  up  to  32  processors 
with  a  maximum  grid  resolution  of  4,370,000  vertices  were  performed  with  KIVA-3/MPI 
(acronym  for  the  present  parallel  version  of  KIVA-3).  A  Smagorinsky  type  subgrid-scale  (SGS) 
model  was  facilitated  during  these  simulations.  Also  an  improved  convection  scheme 
combination  was  used  by  introducing  a  third  option  into  KIVA-3  which  simply  selects  a 
convection  scheme  based  on  central  differencing  for  the  convective  terms  in  the  momentum 
equations  and  Quasi-Second-order  Upwind  (QSOU)  for  the  density  and  energy  equations  (Gel  et 
al.  (1998)).  Up  to  3  %  grid  stretching  was  employed  in  radial  plane  (y-z)  to  capture  the  core 
flow  and  maintain  the  second  order  spatial  accuracy  as  close  as  possible. 

Several  preliminary  large-scale  runs  over  one  million  vertices  were  performed  for  the 
second  and  third  cases  at  different  grid  resolutions.  Table  4.3  shows  the  grid  specifications,  the 
total  number  of  processors  (PEs)  employed,  and  the  average  CPU  time  required  per  each 
timestep  of  the  simulation  for  the  second  case  (i.e.,  turbulent  round  jet).. 


Table  4.3  CPU  seconds  per  timestep  for  the  preliminary  production  run 


Tot.  Vertices 

#PEs 

CPU 

160  X  80x80 

1,089,288 

16 

32.0 

160  X  80x80 

1,089,288 

32 

17.0 

208  X  100  X  100 

2,184,840 

16 

62.0 

288  X  122  X  122 

4,370,000 

32 

130.0 

Single  processor  KIVA-3  runs  at  these  grid  resolutions  that  are  required  for  the 
calculation  of  the  speedup  and  parallel  efficiency  for  these  cases  could  not  be  performed  due  to 
very  long  turn-around  times.  However,  as  seen  from  Table  4.3  for  the  first  two  cases,  when  the 
number  of  processors  employed  were  doubled  for  the  same  size  problem,  the  average  CPU  time 
per  timestep  nearly  reduced  in  half  which  demonstrates  a  nearly  linear  speedup. 
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The  third  problem  was  the  free  swirling  annular  jet  case  studied  by  Smith  et  al.  (1999a). 
An  annular  jet  with  an  inner  and  outer  radius  of  2.7  cm  and  5.3  cm,  respectively,  enters  a  three- 
dimensional  rectangular  domain  from  the  left  wall  with  an  initial  swirl  velocity  profile.  Initial 
conditions  were  based  on  the  measurements  of  Holzapfel  (1996).  Table  4.4  illustrates  the 
speedup  and  parallel  efficiency  achieved  for  the  80  x  80  x  80  grid  configuration  up  to  16 
processors.  Note  that  the  numbers  reported  here  are  based  on  setup  phase  excluded  timing.  Other 
benchmark  results  can  be  found  in  Gel  (1999),  Gel  &  Celik  (1999a,  1999b,  1999c). 

Table  4.4  Speedup  &  Efficiency  for  80  x  80  x  80  (500K)  Grid 


Number  of  PEs 

Speedup 

Efficiency 

1 

1.0 

1.00 

2 

1.7 

0.86 

8 

6.6 

0.82 

16 

11.4 

0.71 

4.3.2  Hardware  Improvements 

As  part  of  an  effort  to  improve  the  local  computational  resources,  a  20  node  DEC  Alpha 
cluster  was  built  within  our  CFD  Laboratory.  Currently  12  nodes  are  active  but  the  other  nodes 
can  be  easily  activated.  Each  compute  node  consists  of  533  MHz  DEC  Alpha  21164  processor 
with  256  MB  RAM,  8.3  GB  harddisk,  two  DEC  Tulip  chipset  based  Kingston  TXlOO  Fast 
Ethernet  cards.  The  head  node  consists  of  the  same  configuration  except  that  the  memory  is  512 
MB  and  an  additional  ethemet  card  is  installed  for  accessing  the  Beowulf  cluster  from  outside 
networks  to  submit  jobs.  The  operating  system  is  the  Red  Hat  Linux  5.2  with  2.0.35  kernel  level. 
Each  system  is  configured  to  boot  standalone  without  the  necessity  of  a  keyboard,  a  mouse  and  a 
graphics  card  to  be  attached.  The  system  configuration  changes  are  performed  through  serial  port 
where  the  I/O  at  the  reboot  is  being  forwarded.  This  configuration  reduced  the  cost  of  each  node 
significantly.  In  order  to  achieve  a  peak  I/O  throughput  of  200  Mbit/s  via  channel  bonding,  the 
nodes  were  connected  to  each  other  through  two  24  port  Fast  Ethemet  switches. 
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Table  4.5  Speedup  and  Efficiency  Comparison  for  80  x  80  x  80  Grid 


No.  of  Processors 

SoeeduD 

Efficiency 

Origin  2000 

1 

1.0 

1.00 

2 

1.7 

0.86 

8 

6.6 

0.82 

16 

11.4 

0.71 

Beowulf  Cluster 

1 

1.0 

1.00 

8 

5.1 

0.64 

10 

6.0 

0.60 

Several  benchmark  runs  were  performed  with  the  annular  swirling  jet  test  case  using 
KIVA-3/MPI.  Table  4.5  shows  the  speedup  and  efficiency  figures  this  case.  As  seen  in  this  table, 
the  performance  of  the  new  Beowulf  cluster  is  almost  as  good  as  SGI  Origin  2000.  Further 
details  about  the  cluster  configuration  and  the  benchmark  runs  can  be  found  in  Gel  &  Celik 
(1999d). 
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Chapter  5:  RANS  Modeling  and  Turbulence  Scales 


As  part  of  the  objectives  of  this  project  the  commonly  used  Reynolds  averaged  Navier- 
Stokes  (RANS)  turbulence  models  are  studied  and  the  results  are  assessed  with  regards  to  their 
performance  for  predicting  mean  flow  and  turbulence  scales.  This  study  has  enhanced  our 
understanding  of  the  deficiencies  of  the  RANS  models  and  provided  insight  as  to  how  they  can 
be  remedied.  It  was  also  important  to  establish  the  order  of  magnitudes  of  the  tmbulence  length 
and  time  scales  in  relation  to  the  engine  geometry  and  speed  to  guide  us  in  determining  the 
minimum  mesh  and  time  resolution  for  LES.  This  section  presents  the  important  findings  of  our 
study  concerning  RANS  models  and  turbulence  scales  as  predicted  by  these  models  for  IC 
engines.  Detailed  information  can  be  found  in  the  papers  referenced. 

5.1  Reynolds  Averaged  Navier-Stokes  Models 

The  literature  on  applications  of  computational  fluid  dynamics  (CFD)  to  internal 
combustion  engines  show  that  the  turbulence  model  of  choice  widely  utilized  are  of  the  two- 
equation  models  based  on  isotropic  eddy  viscosity  concept  (see  Celik  et  al.,  1999,  for  a  review). 
However,  there  is  little  information  in  the  common  literature  with  regards  to  the  validation  and 
performance  comparison  for  various  models  of  this  type.  Our  experience  has  been  that  these 
models  do  not  give  results  consistent  with  each  other  when  applied  to  realistic  in-cylinder  engine 
flows.  It  is  found  necessary  to  perform  a  comparative  validation  study  among  various  models, 
which  might  benefit  the  CFD  community  involved  with  in-cylinder  turbulent  flow  predictions. 

The  current  understanding  of  turbulence  phenomenon  as  it  relates  to  combustion,  engine 
performance,  and  pollutant  emissions  is  far  from  being  satisfactory.  The  combustion  process, 
hence  fuel  efficiency  and  pollutant  formation,  (see  e.g.  Hey  wood,  1987)  is  greatly  influenced  by 
the  turbulence  generated  during  the  intake  stroke  followed  by  the  turbulence  induced  by  the 
geometry  of  cylinder-piston  assembly  during  the  compression-expansion  stroke.  The  control  of 
the  consequences  of  turbulence  can  only  be  possible  with  a  well  tested  and  validated  turbulence 
model  specifically  tuned  for  IC  engine  applications. 

There  is  an  extensive  amount  of  literature  on  RANS  modeling  for  IC  engines  (Bo  et  al., 
1997;  Celik  et  al.,  1992,  1997;  Reitz  and  Rutland,  1991;  Han  et  al.,  1997;  Ramos,  1989; 
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O’Rourke  and  Amsden,  1987;  Bashay  et  al.,  1986;  Gosman,  1985).  Most  of  the  current 
computational  studies  utilize  some  version  of  the  so-called  k-e  turbulence  model.  Because  of  the 
multiple  eddy  structures  and  the  strongly  elliptic  nature  of  the  governing  equations,  low-order 
(i.e.  zero-equation)  models  do  not  work  well.  Some  studies  use  algebraic  stress  models  (ASM), 
but  the  improvements  over  the  k-8  model  and  its  variants  (e.g.  low  Reynolds  number  k-s  model; 
RNG  k-e,  Lien  and  Leschziner,  1994,  Han  et  al.,  1996)  are  marginal  and  there  is  no  consensus  in 
the  literature  as  to  which  model  performs  better.  The  isotropic  eddy  viscosity  concept  and 
gradient  diffusion  model,  where  the  turbulent  stresses  are  related  solely  to  the  local  strain  rates, 
are  central  to  the  k-s  model.  It  is  well  known  (see  e.g.  Reynolds,  1980)  that  in  most  cases 
relevant  to  engine  combustion  chambers  the  above  assumptions  are  not  correct.  Due  to  this  fact, 
the  performance  of  the  k-s  model  has  not  been  fiilly  successful,  and  the  predictions  are  not  as 
accurate  as  it  is  desirable;  in  some  cases  performance  is  very  poor.  The  comparisons  presented 
by  Leschziner  (1991),  Celik  et  al.  (1987),  and  Mongia  et  al.  (1986),  confirm  the  above 
arguments.  The  applications  of  second  moment  closure  models  (i.e.  Reynolds  stress  transport- 
RST  models)  to  IC-engines  are  scarce  in  the  literature  for  the  obvious  reason  that  six  additional 
non-linear  partial  differential  equations  must  be  modeled  and  solved.  RST  models  were  applied 
to  engine-like  geometries  in  a  recent  study  by  Yang  et  al.  (1998).  Although  this  study  shows 
some  promise  it  remains  to  be  validated.  This  study  also  shows  significant  differences  between 
k-e  and  RST  model  results.  Leschziner  (1991)  has  demonstrated  that  while  considerable 
improvement  can  be  gained  from  second  moment-closure,  predictive  performance  is  not 
consistently  satisfactory.  Hence  the  two-equation  turbulence  models  remain  the  modelers  choice 
in  IC  simulations,  and  their  predictive  capability  need  to  be  assessed. 

In  the  present  study,  an  isothermal,  incompressible  flow  within  a  piston-cylinder 
arrangement  (Morse  et  al.,  1979;  Haworth,  1998)  motored  without  compression  at  200  RPM  is 
investigated  (Yavuz  and  Celik,  1999a,b).  The  mean  piston  speed  Vp  was  40  cm/sec  in  these 
calculations.  There  was  no  swirl  imposed.  The  seat  angle  was  60*  with  the  horizontal  axis.  A 
time  step  of  1x10’^  seconds  was  used  for  time  marching.  A  numerical  mesh  of  40,000  vertices 
was  employed  with  a  wedge  angle  of  0.5  degrees  in  the  circumferential  direction.  The  maximum 
cell  size  for  the  wedge  calculations  was  0.5x0.5  mm  in  the  radial  and  axial  direction  inside  the 
cylinder.  The  inlet  pressure  was  atmospheric  and  the  initial  flow  was  at  rest  till  the  piston  started 
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to  move  away  from  TDC  (0°  CA).  Five  models  were  considered,  namely,  (1)  The  standard  k-e 
model,  the  RNG-k-e  model,  a  low  Re  number  k-e  model,  a  modified  k-e  model  for  curvature 
effects  (RNC-k-e;  see  Yavuz  and  Celik,  1999a),  and  a  new  Smagorinsky  based  eddy  viscosity 
(SEV;  see  Yavuz  and  Celik,  1999b)  model.  All  of  these  were  implemented  into  the  KIVA  code. 

The  predicted  streamlines  from  the  intake  flow  case  are  shown  in  Figure  5.1  as  compared 
to  the  measurements  of  Morse  et  al.,  (1979).  The  results  from  all  models  compare  well  with  the 
measurements  in  the  overall  sense.  The  RNG  k-e  model  solutions  deviate  more  from  the 
measurements  in  capturing  the  size  of  the  lower  wall  recirculation  zone.  The  Low-Re  k-e  model 
performs  slightly  better  in  this  regard,  which  is  supported  by  the  fact  that  the  turbulent  Reynolds 
number  is  low  in  such  a  flow. 

In  Figure  5.2  the  axial  velocity  profiles  are  compared  at  various  stations  inside  the 
cylinder.  The  standard  k-e  model  and  the  Low-Re  model  results  are  essentially  the  same  and 
compare  well  with  the  measurements.  The  Low-Re  model  seems  to  improve  the  results  near  the 
wall  (as  expected)  as  well  as  in  the  vicinity  of  the  reattachment  point  at  the  wall. 

Then  the  flow  field  inside  a  typical  piston-cylinder  assembly  (Catania  et  al.,  1995)  was 
simulated  under  motored  conditions  at  600  RPM,  which  implies  a  mean  piston  speed  Vp  =150 
cm/sec.  A  numerical  mesh  of  20,000  vertices  was  used  with  a  wedge  angle  of  0.5  degrees  in  the 
circumferential  direction.  The  maximum  cell  size  for  the  wedge  calculations  was  0.6x0.9  mm  in 
the  radial  and  axial  direction.  As  reported  by  Celik  and  Yavuz  (1997),  the  minimum  integral 
length  scale  occurs  at  the  end  of  the  compression  stroke  in  the  order  of  squish  clearance  height. 
To  investigate  grid  dependency  another  mesh  with  70000  vertices  was  used  for  the  same  case, 
which  provided  approximately  the  same  results  in  the  mean  flow  for  the  standard  k-s  model. 
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Figure  5.2  Profiles  of  axial  velocity  at  90°  CA  (Yavuz  and  Celik,  1999a,  1999b) 


The  results  of  the  standard  k-e  and  RNG  k-e  models  are  compared  the  piston-bowl 
assembly  (Figure  5.3).  While  no  significant  differences  were  observed  in  the  intake  flow  case, 
solutions  for  the  piston-bowl  assembly  are  drastically  different  at  some  crank  angles.  That  again 
raises  the  question,  if  these  models  are  applicable  to  IC-engine  flows  at  all  speeds  and  operating 
conditions?  Even  if  they  both  predict  good  results  for  one  case,  they  could  predict  significantly 
different  results  for  other  cases,  which  immediately  challenges  the  validity  of  the  solution. 


57 


Figure  5.3  Streamlines  for  the  piston-bowl  case  at  90°  CA  (Yavuz  and  Celik,  1999a,  1999b) 

In  general  all  of  these  models  do  a  reasonably  good  job  in  predicting  the  main  features  of 
the  mean  flow  during  the  intake  stroke  when  the  valve  is  not  moving  and  the  piston  speed  is  low. 
Detailed  quantitative  comparisons  revealed  that  the  Low-Re  model  and  the  RNC  k-e  (with 
curvature  correction)  model  seems  to  perform  the  best,  with  an  appreciable  improvement  near 
the  walls  for  the  case  where  the  relatively  low  turbulence  Reynolds  numbers  were  observed.  In 
this  case  the  RNG  k-s  model  performs  worse  than  the  standard  k-e  model  contrary  to  our 
expectation.  Perhaps  a  low  Reynolds  number  version  of  this  model  would  perform  better.  The 
SEV  model  has  the  potential  to  predict  the  physics  of  this  flow  better,  provided  that  the 
characteristic  length  scale  is  calculated  dynamically,  including  all  the  relevant  scales  of  the 
engine,  and  the  distance  from  the  wall.  Such  a  model  is  currently  under  investigation  by  the 
authors.  However,  although  the  mean  flow  predictions  are  reasonably  good,  the  turbulence 
intensity  predictions  do  not  compare  well  with  the  experimental  data,  which  would  show  the  real 
performance  of  a  turbulence  model.  An  improved  version  of  this  model  is  being  developed  as  a 
hybrid  Smagorinsky  turbulence  model,  which  is  used  for  flow  predictions  ranging  from  RANS  to 
LES  automatically  depending  on  the  grid  size.  Further  information  can  be  found  in  Yavuz 
(2000). 

As  for  the  case  of  the  piston-bowl  assembly  with  closed  valves  (i.e.  compression  and 
expansion  stroke)  the  two  models  (standard  k-e  and  RNG  k-e)  produce  drastically  different 
results,  even  though  their  performance  was  very  similar  in  the  case  with  intake  flow. 
Unfortunately,  the  resolution  of  this  unexpected  behavior  requires  more  validation  studies  at 


higher  engine  speeds  and  analysis  of  the  results  from  several  engine  cycles  against  a  reliable  set 
of  experimental  data.  It  seems  that  for  IC  flows  transient  effects  are  very  important  and  the 
turbulence  models  need  to  be  tested  against  a  good  time-dependent  (in  the  mean)  benchmark 
case. 


5.2  Turbulence  Scales 

There  is  extensive  amount  of  information  available  in  the  literature  concerning  turbulence 
intensity,  length  scales  and  turbulence  in  general.  However  this  information  which  come  from 
discrete  experiments  which  are  engine  or  problem  specific,  have  rarely  been  assessed  together 
under  the  realm  of  theory  and  engine  simulations  using  turbulence  models.  We  have  reviewed 
only  the  most  relevant  literature  (see  Celik  and  Yavuz,  1997).  Here  we  present  a  summary  of  our 
findings. 

The  relevant  length  and  time  scales  pertaining  to  both  the  geometry  and  to  turbulence  are 
reviewed  by  Reynolds  (1980)  and  Hey  wood  (1987),  and  computed  by  Shah  and  Markatos 
(1986).  The  length  scales  exhibit  a  wide  spectrum  of  sizes  from  the  size  of  the  combustion 
chamber  down  to  the  Kolmogorov  microscale.  In  between  are  the  thickness  of  the  wall  boundary 
layers,  the  width  of  the  intake  jet,  the  integral  length  scale  of  turbulence,  and  the  Taylor 
microscale  as  relevant  length  scales  which  must  be  considered.  The  important  time  scales  include 
the  integral  time  scale  of  turbulence,  large  eddy  turnover  time,  the  flow-pass  time  (a  measure  of 
time  required  for  a  fast  moving  fluid  eddy  to  travel  the  length  of  the  shear  layer),  and 
development  time  (time  required  for  a  boundary  layer  to  reach  equilibrium).  Since  burning 
occurs  at  the  smallest  scales,  it  is  desirable  that  these  scales  are  either  resolved  or  modeled  (e.g. 
via  sub-grid  scale  models)  to  accurately  predict  the  flame  structure  and  propagation.  The  time 
scale  associated  with  the  smallest  eddies  is  the  Kolmogorov  time  scale  (=(v/e)'^  ,  v  being  the 
kinematic  viscosity,  and  e  the  dissipation  rate  of  turbulent  kinetic  energy,  k). 

Hey  wood  (1987)  presents  a  review  of  the  turbulence  length  and  time  scales  in  IC  engines 
which  is  summarized  in  Table  5.1.  An  important  information  pertaining  to  the  length  scales  in 
engines  is  that  the  maximum  wall  boundary  layer  thickness  on  the  cylinder  head  and  piston 
surfaces  could  be  as  high  as  5  mm  at  moderate  engine  speeds.  The  flow  regimes  inside  the 
boundary  layer  as  well  as  the  main  vortex  motion  inside  the  cylinder  can  be  characterized  by  a 
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Reynolds  number  defined  as  Re  =  VpL/v  ;  where  Vp  is  the  mean  piston  speed,  L  is  the  stroke,  and 
V  is  the  kinematic  viscosity. 


Table  5.1  Estimated  length  and  time  scales  for  an  automotive  size  engine  at  1000  RPM 


Estimated 

Value 

Process,  location 

Comments 

10.0  mm 

During  intake 

scales  with  valve  lift 

ll 

(min) 

2.0-5.0  mm 

end  of  compression 

scales  with  clearance  height, 
varies  little  with  engine  speed 

It 

1.0  mm 

end  of  compression 

varies  little  with  engine  speed 

Ik 

0.01  mm 

end  of  compression 

1.0  ms 

end  of  compression 

decreases  with  engine  speed 

Tt 

0.1  ms 

end  of  compression 

decreases  with  engine  speed 

Notation:  h  =  Integral  length  scale  ;  It  =  Taylor  micro  length-scale;  Ik  =  Kolmogorov 
length  scale;  x\i  =  Integral  time  scale;  tx  =  Taylor  micro  time-scale 

Fraser  et  al.  (1986)  presents  direct  measurements  of  lateral  integral  length  scale  (Figure 
5.4)  using  laser  Doppler  velocimetry  (LDV).  This  length  scale  ranges  between  2  to  3.5  mm  for 
crank  angles  between  320  to  380  degrees  ATDC  with  a  local  minimum  at  the  TDC.  When  the 
same  data  is  normalized  by  the  instantaneous  clearance  height  the  non-dimensional  value 
increases  from  0.1  to  0.2  at  TDC  and  remains  approximately  fixed  after  that.  Lancaster’s 
measurements  agree  with  those  by  various  other  authors  within  a  factor  of  two.  The  previous 
work  reviewed  by  Lancaster  along  with  his  work  seem  to  indicate  that  the  relative  integral  length 
scale  of  turbulence  should  be  function  of  the  engine  speed,  swirl  ratio,  and  compression  ratio.  He 
also  points  out  that  the  longitudinal  and  lateral  length  scales  could  be  significantly  different;  for 
homogeneous,  isotropic  turbulence. 

Recently,  Corcione  et  al.  (1990)  presented  results  from  direct  measurements  of  lateral 
length  scale  inside  the  combustion  chamber  of  a  diesel  engine  during  the  compression  stroke  at 
600  and  1000  RPM.  They  calculate  indirectly  the  lateral  length  scale  using  Taylor’s  hypothesis, 
however  due  to  the  uncertainty  involved  in  this  analysis  we  only  consider  the  direct 
measurements.  Their  correlation  coefficients  are  also  of  the  form  of  Equation  (1).  The  measured 
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lateral  length  scales  were  in  the  order  of  1  to  2.5  mm  and  decreases  during  the  compression 
stroke  (330  to  370  degree  crank  angles).  These  authors  showed  that  the  cylinder  geometry  had  a 
significant  influence  on  turbulence  length  scales. 


Detailed  measurements  of  root  mean  square  (rms)  turbulent  velocity  fiuctuations  using  a 
state-of-the-art  LDV  system  are  reported  by  Arcoumanis  et  al.  (1994)  for  a  1.9  liter  direct 
injection  (DI)  diesel  engine.  The  rms  values  of  the  swirl  velocity  fluctuations  are  very  high 
(about  4Vp,  Vp  being  the  mean  piston  speed)  during  the  intake  due  to  counter-rotating  jets,  they 
decay  to  about  IVp  near  the  BDC.  The  rms  axial  velocity  fluctuations  attain  peak  values  of  about 
4Vp  during  early  induction.  Both  swirl  and  axial  velocity  fluctuations  show  as  much  as  three  fold 
spatial  variation  with  radial  position  at  a  dimensionless  distance  x/S  =  0.1  (here  x  is  the  axial 
distance  and  S  is  the  stroke)  from  the  cylinder  head  indicating  that  the  turbulence  is  strongly  non- 
homogeneous. 


Figure  5.4  compares  some  results  from  our  own  computations  (Celik  and  Yavuz,  1997) 
for  an  axisymmetric  flow  during  the  compression  and  expansion  stroke  of  a  typical  IC  engine. 
The  results  indicate  that  the  standard  k-e  model  does  seem  to  give  the  right  orders  of  magnitude 
and  even  the  right  trends  for  the  turbulence  velocity  and  length  scales.  Needless  to  say  that  this 
model  assumes  local  isotropy,  hence  it  can  not  predict  the  usually  observed  non-isotropic 
distribution  of  turbulence,  especially  during  the  intake  and  exhaust  strokes.  The  flow  during 
compression  and  expansion  strokes  is  known  to  be  nearly  isotropic. 
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Figure  5.4  Comparison  of  calculated  (Celik  and  Yavuz,  1997)  dimensionless  integral  length  scale 

with  experiment  (Fraser  et  al.,  1986) 
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Figure  5.5  Integral  and  Kolmogorov  length  scales  inside  an  IC-Engine  (Celik  and  Yavuz,  1997) 

with  and  without  combustion. 


The  computational  results  further  indicate  (Figure  5.5.)  that  the  turbulence  intensity  is 
very  sensitive,  but  the  integral  length  scale  is  not  as  much  sensitive  to  combustion.  Interestingly, 
the  Taylor  and  Kolmogorov  length  scales  are  sensitive  to  combustion,  and  this  influence 
manifests  itself  with  a  considerable  lag  with  respect  to  the  increase  in  temperature. 

The  minimum  length  scales  (integral,  Taylor,  and  Kolmogorov)  are  predicted  near  the 
TDC,  and  are  1  to  5  mm,  0.05  mm,  and  0.01  mm,  respectively,  which  are  in  good  agreement 
with  experimental  observations.  In  this  regard  the  k-s  can  be  used  as  a  guide  to  determine 
computational  strategies  for  fine  grid  LES  simulations.  It  seems  that  it  would  be  logical  to  isolate 
the  zone  near  top  dead  center  and  use  many  more  spatial  grid  nodes  in  this  zone  as  compared  to 
other  zones  with  accurate  interpolation  and/or  extrapolation  to  continue  the  computations 
elsewhere. 
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Chapter  6:  Development  and  Application  of  A  Subgrid-Scale  PDF  Combustion 
Model 

6.1  Background 

In  this  chapter  the  development  and  application  of  a  subgrid  scale  combustion  model  to 
be  used  for  Large  Eddy  Simulation  studies  is  presented.  The  theoretical  basis  of  the  model  is 
described  first.  Then,  the  model  capabilities  are  demonstrated  by  comparing  the  combustion 
predictions  to  some  available  experimental  benchmark  cases.  Because  most  of  the  experimental 
measmements  are  available  as  time  averaged  or  as  Favre  averaged  quantities  direct  validation  of 
instantaneous  computations  is  not  possible.  However,  predictions  of  instantaneous  temperature 
in  a  two-dimensional  model  will  be  presented  to  demonstrate  the  capabilities  of  the  proposed 
model.  Before  going  into  the  explanation  of  the  model  development,  we  will  present  first  the 
development  of  some  auxiliary  studies  and  computations  that  were  initially  conducted  to 
investigate  the  capabilities  of  Kiva-3  code  to  be  used  in  Large  Eddy  Simulation  studies. 

6.2  Large  Eddy  Simulation  of  a  Round  Turbulent  Free  Jet 

This  investigation  is  aimed  at  research  in  the  large  eddy  simulation  (LES)  of  a  passive 
scalar  mixing  problem  using  Kiva-3  code.  The  preliminary  simulation  performed  here  is  for  a 
free  round  turbulent  jet  of  methane.  The  jet  is  at  Reynolds  number.  Re  =  2000  and  it  is 
discharging  into  stagnant  air.  This  small  Re  is  chosen  in  order  to  be  able  to  resolve  the  basic 
features  of  free  jet  turbulence  using  the  available  workstation  computer  power.  The  computations 
are  started  from  the  near  orifice  field  with  uniform  initial  tmbulence  intensity  and  length  scales. 
No  subgrid  eddy  viscosity  models  has  been  used  as  the  subgrid  dissipation  is  achieved  through 
numerical  truncation  errors.  The  convection  terms  are  discretized  by  a  hybrid  scheme  with  10% 
upwinding  (i.e.  first  order)  and  90%  central  differencing.  The  second  order  (Crank  Nielson) 
temporal  discretization  scheme  is  used  for  the  pressure  gradient  and  diffusion  terms.  The 
convection  terms  are  discretized  by  a  first  order  accurate  scheme  in  time. 

The  all-Mach  number  formulation  employed  in  Kiva-3  allows  the  use  of  this 
compressible  flow  code  to  simulate  the  weakly  compressible  conditions  by  employing  the  proper 


63 


pressure  gradient  scaling.  The  time  step  is  limited  by  courant  stability  requirements  and  it  is  in 
the  order  of  1 0"^  seconds.  The  computational  domain  in  radial  and  azimuthal  directions  is  five 
orifice  diameter  (d)  while  in  axial  direction  it  is  20d.  Outlet  flow  boundary  conditions  with 
allowance  for  flow  reversal  are  assigned  to  the  surrounding  of  the  jet.  This  choice  for  the 
boundary  conditions  is  based  on  the  code  limitations  rather  than  representing  the  actual  situation. 
The  initial  instability  was  introduced  into  the  computations  by  perturbing  the  axial  velocity 
randomly  by  20%  to  simulate  fluctuations  due  to  turbulence  in  the  nozzle.  In  the  latter  stages  of 
the  computations  this  artificial  forcing  is  removed  to  see  whether  the  computed  flow  instabilities 
could  be  self-sustained. 

The  computations  are  continued  several  flow-through  times  which  is  defined  as  tf  =  L/Uo, 
where  L  is  the  length  of  the  domain  and  Uo  is  the  mean  flow  bulk  velocity.  Figure  6.1  shows  the 
transition  from  the  inviscid  potential  core  to  the  chaotic  turbulent  structure  to  occur  at 
approximately  3-4d  downstream  of  the  jet  inlet.  At  the  early  stages  in  the  computations  (t<  tf) 
this  location  is  strongly  affected  by  the  development  of  the  large-scale  structures  in  downstream 
locations  (Grinsten,  et  al.,  1987).  It  is  known  that  the  near-field  of  the  free  jet  is  also  affected  by 
the  large-scale  structures  which  form  in  the  annular  region  between  the  jet  and  the  surroundings 
as  a  result  of  the  shear  layer  instability  (Gladnick,  et  al.,  1990).  This  feedback  process  from 
organized  large  eddies  affects  the  shear  layer  ahead  of  the  nozzle  for  times  less  than  the  flow¬ 
through  time  and  continues  to  do  so  after  several  flow-through  times.  This  feedback  is  also 
reflected  in  the  decrease  of  the  time  step,  which  is  computed  dynamically  from  Courant  number 
limitations  during  computations.  In  the  computations,  for  reason  of  numerical  accuracy  (and 
stability)  the  Courant  number  is  set  to  0.1  of  its  maximum  value.  The  details  of  these  simulations 
can  be  found  in  Amin  et  al.,  (1998a). 

Figures  6.1  and  6.2  show  the  sequential  unsteady  events  captured  with  Kiva-3 
computations.  The  captured  flow  events  such  as  transition  from  laminar  to  turbulent  flow,  vortex 
pairing  and  breakdown,  and  shear  layer  instability  are  in  good  agreement  with  Frankel  et  al.’s 
(1993)  computations. 

The  round  jet  simulations  show  that  the  difference  schemes  might  have  an  important 
influence  over  the  ability  of  the  computed  flow  instabilities  to  persist.  The  use  of  the  flux 
limiting  scheme  known  as  Quasi  Second  Order  Upwinding  (QSOU),  although  it  approaches 
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second  order  accuracy,  may  lead  to  disappearance  of  the  diffusive  wiggles.  Thus  the  flow 
instabilities  will  not  be  able  to  persist.  Higher  grid  resolution  may  require  introduction  of 
artificial  eddy  viscosity  or  changes  in  the  properties  of  the  advection  scheme  in  order  to  maintain 
computational  stability. 


6.3  Probability  Density  Function  Model  Development 

In  mixing  limited  combustion  the  rate  of  reaction  is  determined  by  the  mixing  frequency 
between  the  fuel  and  oxidizer  eddies.  The  use  of  conventional  RANS  techniques  fails  to  directly 
capture  the  dynamics  of  mixing-controlling  eddies.  With  the  use  of  Large  Eddy  Simulation 
(LES)  technique  we  are  able  to  describe  accurately  the  large  scale  mixing  process  while  the 
small  scales  are  left  to  be  modeled.  Extension  of  LES  techniques  devised  for  isothermal 
turbulence  to  combustion  problems  is  not  trivial.  This  is  because  of  the  fact  that  chemical 
reactions  occur  within  thin  zones  much  smaller  than  the  grid  scale.  Therefore  knowledge  of  the 
reactants  distribution  in  the  grid  cell  is  an  important  perquisite.  One  strategy  to  account  for  the 
strong  non-linear  reactant  distribution  within  the  grid  cell  is  to  employ  a  probability  density 
function  (pdf)  for  the  mixing  process. 

In  this  section  we  will  present  the  development  of  a  pdf  model  for  combustion  of  direct 
injection  diesel  engines  which  are  characterized  by  mixing  limited  reactions. 

6.3.1  Filtered  mixture  fraction  and  concentration  fluctuations 


The  transport  of  a  filtered  mixture  fraction,  field  can  be  described  by  the  following 
equation. 


dt  3xj  5Xj  <T|  5Xj 


(6.1) 


Where  v  and  x  are  the  velocity  and  position  in  a  space  coordinate  direction  (a=  1,2,3),  pt 
is  the  turbulent  eddy  viscosity  and  is  the  turbulent  Schmidt  number  for  To  be  able  to 
describe  the  statistics  of  the  subgrid-scale  scalar  fluctuations  a  model  for  the  Probability  Density 
Fimction  (pdf)  of  the  scalar  has  to  be  found.  One  of  the  most  commonly  used  approaches  is  to 
assume  the  shape  of  the  pdf  and  then  determine  the  first  two  statistical  moments  of  the 
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distribution  from  modeled  transport  equations  (e.g.  Bilger,  1976,  Amin  and  Celik  (1999b).  The 
transport  equation  of  the  second  statistical  moment  can  determined  from  a  modeled  SGS 
concentration  fluctuations  can  be  written  as  (Spalding,  1971); 


d{p%)  ,  g(/^Vjg) 
5t  5xj 


_a_ 

dx: 


CT^dx.^ 


-  Cgi//, 


5x  • 


n2 


CgzPgf/k 


(6.2) 


The  numerical  values  of  the  constants  in  the  above  equations  are  listed  in  Table  6.1.  In 
the  subsequent  analysis  we  follow  a  one-equation  model  approach  in  which  k  is  determined  form 
a  modeled  equation  and  e  is  determined  form  a  prescribed  grid  dependent  length  scale. 


Table  6.1  Values  of  constants  in  the  k-e-g  model  of  turbulence. 


Cgi 

Cg2 

^g 

bo 

1.84 

0.7 

0.7 

6.3.2  The  large  eddy  probability  density  function 

The  large  eddy  probability  density  function  (pdf)  gives  a  distribution  of  mixing 
frequencies  and  therefore  we  are  not  restricted  to  a  single  mixing  frequency  as  in  the  case  of 
conventional  time  averaging  techniques.  In  the  past  many  experimental  studies  have  been 
conducted  to  determine  the  most  appropriate  pdf  shape  to  represent  the  grid  cell  mixing  process. 
Various  pdf  shapes  that  can  be  found  in  typical  flows  are  discussed  by  (Bilger,  1989).  The  most 
commonly  used  pdf  forms  are  the  Gaussian  (Lockwood  and  Naguib,  1975),  the  Beta  function 
(Jancika  and  Kollman,  1978)  and  the  Dirac-Delta  functions  (Khalil,  et  al,  1975).  The  use  of  the 
Clipped  Gaussian  pdf  is  justified  on  the  bases  of  experimental  observations  and  Direct 
Numerical  Simulation  (DNS)  data.  These  studies  indicate  that  the  passive  scalar  pdf  at  the  final 
stages  of  mixing  is  relaxed  to  a  Gaussian  distribution.  However,  the  initial  stages  of  mixing, 
which  can  be  approximated  by  a  Double  Delta  function,  cannot  be  well  represented  by  assuming 
a  fixed  Gaussian  form  for  the  pdf  throughout  the  mixing  process.  Beta  function  seems  to  be  able 
to  better  represent  the  different  mixing  stages  including  the  bimodal  distribution.  In  addition,  it 
has  the  potential  for  extension  to  multiple  scalar  mixing  (Girimaji,  1991).  It  has  been  also  used 
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in  many  previous  LES  studies  (e.g.  Cook  and  Riley,  1994).  The  p  (^)  is  assumed  here  to  follow 
a  Beta  function  distribution  which  is  given  by; 


(6.3) 


The  exponents  a  and  p  are  related  to  the  Favre  filtered  mixture  fraction  and  the  subgrid- 
scale  variance  of  the  turbulent  part  of  the  composite  pdf. 


Thus,  the  filtered  density  can  be  determined  from; 

P(x)  =  ipfeW«.x)d^  (6.6) 

0 


6.3.3  Large  eddy  species  mass  fraction  and  reaction  rate 

Two  of  the  main  objectives  of  LES  of  reacting  flows  are  the  determination  of  Large  Eddy 
Mass  Fraction  (LEMF)  and  Large  Eddy  Reaction  Rate  (LERR).  The  LEMF  can  be  determined 
from 

Y,=4i;Y,fe)Pfe)di  (6.7) 

p 

Similarly  the  LERR  can  be  obtained  from 

»i=4);w,fe)pfe>i^  (6.8) 

P 

In  the  above  equations  the  relationships  Yi  (^)  and  Wj  (^)  are  thermo-chemical  state- 
relationships  obtained  from  the  fuel/oxidizer  reaction  chemistry. 

Because  of  the  inherent  difficulty  in  obtaining  the  thermo-chemical  state-relationships  as 
a  function  of  local  mixing  intensity  and  other  external  factors  (e.g.  varying  pressure  conditions. 
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radiative  heat  loss,  turbulence/radiation  interaction)  integration  of  expressions  6.7  and  6.8  could 
represent  a  demanding  computational  effort.  Therefore  a  more  simplified,  yet  accurate  approach 
seems  to  be  in  favor  for  the  3-D  modeling  of  engines  combustion  chambers.  In  the  next  section 
we  present  the  development  and  application  of  a  simplified  model  to  determine  the  large  eddy 
reaction  rate  which  is  originated  from  the  commonly  used  Eddy  Breakup  Model  (EBU).  This 
model  is  widely  used  by  the  combustion  modeling  community  in  general  and  the  Kiva  code  users 
in  particular. 

6.4  A  Modified  Eddy  Breakup  Model  For  Turbulent  Combustion  Modeling  in 
Diesel  Engines 

In  this  section  we  present  an  application  of  a  modified  form  of  the  Eddy  Breakup  (EBU) 
turbulent  combustion  model  originally  developed  by  Brian  Spalding  of  the  Imperial  College 
more  than  two  decades  ago.  The  later  modifications  allows  the  EBU  model  coefficient  to  be 
dynamically  linked  to  the  flame  surface.  This  link  is  achieved  through  the  use  of  a  Conditional 
Probability  Density  Function  (CPDF).  The  application  of  the  model  is  not  computationally 
intensive  and  can  be  easily  incorporated  in  multidimensional  CFD  codes  like  Kiva-3.  This  makes 
it  suitable  to  be  used  as  a  subgrid-scale  turbulence/chemistry  interaction  model  for  large  eddy 
simulation  of  turbulent  combustion  in  diesel  engines.  The  application  of  the  modified  EBU 
model  to  the  case  a  D.  I.  diesel  engine  using  KIVA  III  has  shown  the  importance  of  this 
modification  in  modeling  of  chemistry /turbulence  interaction  for  I.C.  engines  (Amin  et  al., 
1998b). 

In  the  present  analysis  we  will  follow  the  approach  of  Bilger  and  coworkers  in  modifying 
the  EBU  coefficient.  Consider  the  case  of  a  single  step  overall  heat  releasing  reaction  of  the 
form; 

VpF  +  VoO^VpP  (6.9) 

Where  F,  0,  and  P  denote  fuel,  oxidizer  and  product,  respectively.  The  coefficient  Vj  is 
the  stoichiometric  coefficient  for  species  i.  The  stoichiometric  oxidizer/fuel  ratio,  r,  can  be  then 
computed  fi'om  the  following  relation; 
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In  the  above  expression  W-,  is  the  molecular  weight  of  species  i. 

The  average  rate  of  diffusion  flame  burning  in  the  EBU  model  is  expressed  as; 

^  =  (6.11) 

The  over-bars  indicate  ensemble  averaging,  and  tilde  indicates  a  density-weighted  filtered 
reaction  rates.  The  symbols  p,  k,  and  8,  are  the  filtered  mixture-density,  the  kinetic  energy  of 
turbulence,  and  its  dissipation  rate  respectively.  In  the  above  formulation  Ymin  stands  for  the 
minimum  of  the  reactant  concentrations; 

Y.,.=min(Y„Yo/r)  (6.12) 

In  expression  6.11  A  is  the  model  coefficient  and  is  taken  here  to  be  5.0.  It  should  be 
noted  however  that  a  wide  range  of  values  exists  for  this  coefficient  (typically  in  the  range  1 .0  to 
7.0).  This  can  be  helpful  sometimes  to  time  the  computed  results  to  that  observed 
experimentally.  However,  in  the  present  investigation  we  are  following  Bilger’s  approach  to 
determine  the  value  of  this  coefficient  as  a  function  of  the  mixing  characteristic  in  the  engine. 
We  define  the  mixture  fraction  (4)  which  is  related  to  the  stoichiometry  of  the  reaction  by; 

^  =  (6.13) 

r  +  l 

In  the  detailed  analysis  shown  by  Bilger  (1976a&b)  the  mean  reaction  rate  can  be  written 
as; 


^  =  -5P.Yb&P;«.)  (614) 


In  the  above  expression 


Where  ^st  is  the  stoichiometric  mixture  fraction,  Yf,i  is  the  fuel  mass  fraction  in  the  fuel 
stream  (in  our  case  it  is  equal  to  unity),  Yo,2  is  the  oxidizer  mass  fraction  in  the  air  stream  (= 


0.232  for  O2  in  air),  and  pst  is  the  mixture  density  of  the  stoichiometric  composition.  The 
quantity  Xst  is  the  mean  scalar  dissipation  rate  evaluated  at  the  flame  sheet  and  is  equal  to  (e.g. 
Bilger,  1976a&b); 

(6.16) 


Where  peff  is  the  effective  viscosity,  Prg  is  the  Prandtl  number  of  the  concentration 
fluctuations  («  0.7)  and  is  the  spatial  gradient  of  the  filtered  mixture  fraction.  The  mean 
scalar  dissipation  rate  is  an  important  quantity  that  is  related  to  the  degree  of  species  mixing  in 
the  flowfield  and  determines  the  status  of  combustion  completeness  in  the  combustor  (Peters, 
1984).  The  final  quantity  in  the  previous  expression  for  mean  reaction  rate  (P^),  is  the 
probability  density  function  (pdf)  for  the  mixture  fraction  The  assumed  shape  of  the  pdf  has 
been  defined  previously  in  expression  6.3-6.5. 


In  order  to  evaluate  the  filtered  square  of  the  concentration  fluctuations,  g;  A  transport 
equation  can  be  solved  for  this  quantity  (Spalding,  1971),  however  we  chose  here  to  simplify  the 
analysis  by  assuming  that  the  local  value  of  filtered  g  can  be  evaluated  by  the  equilibrium 
between  production  and  dissipation  of  the  concentration  fluctuations.  Thus  the  local  value  of  g 
can  be  computed  from  the  following  algebraic  expression; 


(6.17) 


Where  Cg  is  a  constant  of  order  0.1.  Computing  g  this  way  will  give  a  good 
approximation  to  the  level  of  concentration  fluctuations  without  having  to  solve  a  separate 
transport  equation.  This  simplification  reduces  the  computational  cost  especially  in  demanding 
applications  of  LES  modeling  of  turbulent  reacting  flows.  To  evaluate  k  and  e;  a  transport 
equation  can  be  solved  for  k  while  8  can  be  computed  using  a  prescribed  length  scale.  The  large 
eddy  reaction  rates  computed  by  the  above  expressions  were  coded  and  linked  to  the  Kiva-3 
code. 

The  coefficient  A  is  computed  dynamically  during  the  solution  procedure  by  equating  the 
large  eddy  rate  expression  6.1 1  and  6.14.  The  dynamical  value  of  this  coefficient  is  termed  here 
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as  Ad.  The  large  eddy  reaction  rates  computed  from  Expression  6.3  and  6.6  will  be  called  here 
cb  pdf  and  <i)EBu  respectively. 


Ad=A- 


(6.18) 


Thus  the  coefficient  Ad  is  equal  to; 

2  % 


(6.19) 


Where  the  ratio  between  the  large  eddy  mass  fractions  can  be  expressed  as  (Bilger; 


1989); 


=  J.(zs.)Vf 


(6.20) 


Here  we  make  an  assumption  based  on  the  experimental  data  given  by  Mudford  &  Bilger 
(1987),  for  the  function  Ji  as; 


J,(Z,)  =  0.45exp(|-zJ 
And  Zst  is  defined  as; 


(6.21) 


(6.22) 


After  some  manipulation  this  leads  to  the  following  expression  for  the  dynamic 
coefficient  Ad; 


A.  P.fe) 

P  J,(zj 


(6.23) 


The  above  value  of  Ad  is  used  in  the  present  subgrid-scale  modeling  of  turbulent 
combustion  in  a  diesel  engine  configuration. 
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6.5  Validation  of  the  PDF  model  in  RANS  jet  mixing  studies 

Before  we  use  the  developed  pdf  model  in  LES  of  reacting  flow  in  diesel  engine 
geometry  we  conducted  some  validation  numerical  experiments  and  compared  the  computations 
with  some  published  measurements  of  the  pdf  in  free  jet  mixing  of  propane  and  air.  This  test  case 
was  studied  experimentally  at  SANDIA  National  Lab.  The  key  parameters  are  listed  in  Table 
6.2.  The  comparison  between  the  predictions  of  the  assumed  pdf  shape  is  shown  in  Figs.  6.4-6. 7. 
The  variation  in  calculated  pdf  along  the  jet  centerline  (i.e.  Y/D  =  0.0)  in  the  initial  mixing 
region  (X/D  =3.08)  is  shown  in  Fig.  6.4.  The  measured  pdf  is  more  biased  towards  the  fuel  side 
and  the  computed  pdf  at  this  location  is  flatter  and  wider  than  measurements.  This  indicates  that 
some  degree  of  mixing  is  achieved  at  this  location  (i.e.  a  faster  jet  decay  is  predicted).  Further 
downstream,  (at  X/D  =  24.78)  the  experimental  pdf  shows  a  clear  Gaussian-like  distribution  and 
slightly  over-predicting  the  propane/air  mixing  (Fig.  6.5).  Figures  6.6  and  6.7  compare  the 
evolution  of  the  pdf  at  two  different  axial  locations  (X/D=15,  X/D=50)  and  the  same  radial 
location  of  Y/D  =  0.55.  Again  the  predicted  pdf  seems  to  overestimate  the  mixing  at  small  axial 
distances  while  further  downstream  the  comparison  between  the  predicted  and  measured  pdf  is 
satisfactory.  In  the  mixing  region  the  distribution  is  bimodal  and  consists  of  a  contribution  from 
both  the  unmixed  air  and  mixed  propane  and  air.  Outside  the  mixing  region  the  assumed  pdf 
could  not  simulate  the  spiked  natoe  of  mixing  profile  and  severely  underestimates  the 
fluctuations  level  (pdf  width). 


Table  6.2  Test  section  dimensions  and  inlet  conditions  for  propane/air 
jet  mixing  case  of  SANDIA. 


Test  Section 

.30  m  X  .30  m 

Jet  Tube  Exit  Diameter 

0.052  m  (I.D.),  0.09  m  (O.D.) 

Jet  velocity  (Bulk) 

53  m/s 

Co-flow  Air  Velocity 

9.2  m/s 

Reynolds  Number 
(based  on  jet  exit  diameter) 

68,000 

Co-flow  Air  Turbulence 

0.4% 

Axial  Pressure  Gradient 

6  Pa/m 
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6.6  Application  of  the  developed  pdf  model  in  combustion  simulation  in  diesel 
engines 

In  this  section  we  will  present  an  application  of  the  developed  pdf  in  LES  of  the  diesel 
engine  case  of  Aoyagi  (1980)  et  ah,  the  simulated  cylinder  and  piston  configuration  listed  in 
Table  6.3.  A  two-dimensional  fine  grid  computational  mesh  composed  of  80  nodes  in  x  direction 
and  250  nodes  in  z  direction  is  used.  The  advection  scheme  accuracy  and  time  accuracy  is 
similar  to  the  ones  used  for  the  free  jet  simulation  in  section  6.2.  The  k-£  sgs  turbulence  model 
was  used  in  this  simulation. 

To  start  with  suitable  initial  condition  for  the  engine  flow  simulation,  a  complete 
isothermal  (motored)  engine  cycle  have  been  simulated  first 

The  computations  show  the  development  of  turbulence  in  the  compression  and  expansion 
strokes  (Figures  6.8-6. 9).  These  figures  indicate  the  chaotic  nature  of  the  turbulent  flow  with 
various  eddy  sizes  distributed  non-homogeneously  throughout  the  cylinder  and  show  great 
variations  with  time.  The  turbulence  generation  did  not  die  out  as  a  result  of  compression  near 
TDC.  As  the  piston  moves  upward  the  maximum  length  scale  reduced  from  1.5  mm  at  -44  deg. 
ATDC  to  0.1  mm  at  -13  deg.  ATDC.  The  velocity  vector  contours  show  that  eddies  of  different 
recirculation  strength  is  found  near  the  edge  and  comer  of  the  piston  walls.  During  expansion 
the  strength  of  the  recirculating  eddies  did  not  die-out  and  disappear.  In  fact,  the  velocity 
magnitudes  show  extensive  localized  recirculation  in  different  parts  of  the  computational 
domain.  During  the  expansion  stroke  the  turbulence  intensity  is  sustained  and  vorticial  strength 
increases.  At  the  BDC  position  (i.e.  180  deg.)  the  contours  of  the  subgrid  kinetic  energy  of 
turbulence  shows  the  different  eddies  swept  by  the  expansion  down  to  the  BDC.  The  velocity 
contours  at  315  deg,  which  is  just  before  the  onset  of  injection  starts  indicated  that  the  flow  is 
quite  turbulent  and  a  random  initial  flow  field  could  be  successfiilly  established.  Figure  6.10 
shows  the  instantaneous  radial  velocity  at  a  point  near  the  cylinder  head.  The  temperature 
variation  at  the  same  point  is  depicted  in  Figure  6.1 1.  Near  the  TDC  the  mean  Radial  velocity 
increased  but  the  fluctuations  did  not  die  out.  In  fact  the  turbulent  fluctuations  intensity  seemed 
to  be  sustained  after  compression  even  without  the  presence  of  the  initial  turbulence  introduced 
by  the  valve  motion  in  the  conventional  engine  geometry.  The  increase  in  mean  temperature  near 
the  onset  of  ignition  is  observed  also  in  these  figures.  Figure  6.10  shows  large-scale  variations  in 
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temperature  as  an  indication  of  unsteady  flame  propagation.  At  about  355°  CA  the  flame  edge 
sweeps  through  this  particular  point,  then  at  about  15°  CA  the  main  flame  front  arrives.  The 
increase  in  velocity  fluctuations  as  a  result  of  combustion  is  noteworthy. 

Table  6.3  Specification  of  the  simulated  engine,  Aoyagi  et  al.  (1980) 


Bore  (mm) 

95 

Stroke  (mm) 

110 

Compression  ratio 

14.6 

Squish  (mm) 

0.9 

Diameter  of  piston  bowl  (mm) 

72 

Height  of  piston  bowl  (mm) 

11.5 

Fuel 

Gas  oil 

Fuel  density  (at  15  de.  C.) 

0.849 

Fuel  cetane  number 

57 

Number  of  injector  nozzle  holes 

4 

Hole  diameter  (mm) 

2 

Injection  angle  (deg.) 

150 

Nozzle  opening  pressure  (MPa) 

17 

Injection  starts  ( °ATDC) 

-15 

Fueling  rate  (mm3/cycle/cylinder) 

39 

RPM 

1250 

Swirl  Ratio 

0 
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Velocity  Vectors  and  Passive  Scalar  Mass  Fraction  Contours 
Round  Turbulent  Free  Jet  After  6  Times  The  Flow-Through  Time 
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Figure  6.1  Intantaneous  contours  of  passive  scalar  mass  fraction  in  a  round  turbulent  free  jet. 
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Figure  6.2  Initial  development  of  the  low  Mach  number  free  jet;  0  <  ti  <  t2  <  .. 


Figure  6.3  Intermediate  development  of  the  low  Mach  number  free  jet  and  initial  development  of 

the  jet  asymmetry;  0  <  ti  <  . . .  <  t?  <  tg  <t9 
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Figure  6.4  Measured  and  predicted  mixture 
fraction  pdf  in  the  jet  (X/D=3.08,  Y/D=0.0). 


Figure  6.5  Measured  and  predicted  mixture 
fraction  pdf  in  the  jet  (X/D=24.78,  Y/D-0.0). 


Figure  6.6  Measured  and  predicted  mixture 
fraction  pdf  in  the  jet  (X/D=l  5,  Y/D=0.55). 


Figure  6.7  Measured  and  predicted  mixture 
fraction  pdf  in  the  jet  (X/D=50,  Y/D=0.55). 
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Figure  6.8  Velocity  magnitude  contours;  (a)  61°  CA,  (b)  183°  CA 
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Figure  6.9  Turbulent  kinetic  energy  contours  near  BDC  (183°  CA) 
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Temperature  (K)  §  Intantaneous  Radial  Velocity  (cm/sec) 


6.10  Instantaneous  radial  velocity  with  Crank  Angle  in  combustion  simulation  at  point 
(x,  y,  z’;  29,  0, 1 17 ),  5.4  mm  below  the  cylinder  head 


Figure  6.1 1  Instantaneous  temperature  with  Crank  Angle  in  combustion  simulation,  at  same 

point  as  in  Figure  6.10 
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Chapter  7  Development  of  Turbulent  Flow  Measurement  Technique 

Several  different  non-intrusive  whole  field  velocimetry  techniques  are  currently  under 
development  which  provide  velocity  data  in  a  plane,  which  can  thus  reduce  the  time  required  to 
map  out  a  complex  flow  field.  This  project  is  exploring  the  accuracy  of  Doppler  Global 
Velocimetry  (DGV),  a  nonintrusive,  planar  imaging,  Doppler-based  velocimetry  technique,  as 
well  as  the  accuracy  of  related  Point  Doppler  Velocimetiy  (PDV).  Both  of  these  techniques  use 
an  iodine  vapor  cell  to  determine  the  Doppler  shift,  and  hence  the  velocity,  of  small  seed 
particles  in  a  flow  field,  as  these  particles  pass  through  a  two-dimensional  sheet  of  laser  light. 
The  same  portion  of  the  light  sheet  is  viewed  through  a  beamsplitter,  either  by  a  pair  of  video 
cameras  (for  DGV),  or  a  pair  of  photodetectors  (for  PDV),  with  the  iodine  cell  placed  in  the 
optical  path  of  one  of  the  cameras  or  photodetectors.  Laser  wavelength  and  cell  absorption  band 
are  matched  such  that  flow  velocities  of  interest  yield  Doppler  shifted  frequencies  in  the  linear 
portion  of  the  absorption  band  of  the  cell.  As  a  result,  the  ratio  of  the  light  intensities  seen  by  the 
two  detectors  at  a  point  in  the  flow  yields  a  signal  that  is  proportional  to  the  particle  velocity. 

A  two-channel  non-scanned  point  PDV  system  has  been  developed  (Kuhlman,  et  al., 
1997),  along  with  a  two-channel  planar  imaging  DGV  system  (Naylor  and  Kuhlman,  1998, 
1999).  The  accuracy  limits  of  both  systems  are  being  systematically  explored,  through  a  series  of 
measurements  in  relatively  simple,  unheated  flows  such  as  fully-developed  turbulent  pipe  flow,  a 
circular  jet,  and  the  flow  over  an  airfoil.  A  rotating  wheel  is  also  being  used  as  a  velocity 
standard.  While  current  DGV  systems  lack  the  accuracy  or  resolution  of  conventional  Laser 
Velocimetry  or  Particle  Image  Velocimetry  (accuracy  of  mean  velocity  typically  on  the  order  of 
5%  of  full  scale,  versus  1%  for  LV),  DGV  has  proven  in  a  short  time  to  be  a  very  flexible  whole- 
field  velocimetry  technique. 

7.1  Apparatus  and  Procedure 

The  present  point  PDV  system  closely  follows  the  basic  DGV  system  that  was  originally 
developed  by  Meyers  et  al.  (1991),  except  that  photodiodes  have  been  used,  along  with  front 
lenses  and  pinholes,  to  collect  scattered  light  from  a  single  point  in  a  seeded  flow,  Kuhlman,  et 
al.  (1997)  and  James  (1997)  have  described  the  PDV  system  in  detail.  A  reference  iodine  cell 
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has  been  used  to  compensate  for  laser  frequency  drift.  Calibration  of  the  iodine  cells  has  been 
accomplished  using  a  continuous  scan  of  the  mode  structure  of  the  cw  Argon  ion  laser,  by 
mechanically  altering  the  tilt  of  the  etalon  through  about  10-20  mode  hops  (James,  1997). 

DGV  system  hardware  and  have  been  described  in  Naylor  and  Kuhlman  (1998,  1999), 
and  in  more  detail  in  Naylor  (1998).  Most  of  this  hardware  is  similar  to  that  used  for  the  PDV 
system.  Eight-bit  Hitachi  KP-Ml  CCD  cameras  and  a  Matrox  Genesis  frame  grabber  have  been 
used  to  capture  images  for  the  DGV  system.  All  four  cameras  are  read  simultaneously  using 
horizontal  and  vertical  sync  pulses  from  the  Genesis  board.  Nikon  35- 135mm,  G. 5-4.5  zoom 
lenses  have  been  used  with  the  CCD  cameras,  because  of  their  versatility  in  imaging  different 
sized  areas  over  a  wide  range  of  distances.  Polarizing  filters  have  been  placed  in  front  of  the 
beam  splitters  to  minimize  effects  due  to  polarization  sensitivity  of  the  beam  splitters.  Image 
processing  software  has  been  developed,  as  described  by  Naylor  (1998),  which  closely  follows 
the  image  processing  methods  developed  at  NASA  Langley  by  Meyers  (1992,  1996).  Flow 
seeding  for  DGV  and  PDV  measurements  is  provided  by  a  commercial  fog  machine,  which  feeds 
a  large  plenum  to  damp  out  pulsations  in  smoke  output. 

Most  recent  results  include  2-component  PDV  data  for  the  flow  over  a  nominally  11.8 
inch  chord  NACA  0012  airfoil  model  at  zero  degrees  angle  of  attack,  as  described  in  Kuhlman 
and  Webb  (1999).  The  airfoil  has  been  mounted  at  the  exit  of  a  2.25x3.25  inch  flow  channel,  fed 
by  a  blower  capable  of  exit  velocities  of  40-45  m/sec.  The  present  results  have  been  obtained  at 
an  exit  velocity  of  21  m/sec,  due  to  limitations  on  removing  the  smoke-seeded  flow  from  the 
laboratory.  Several  series  of  2-component  PDV  data  and  single  component  hot  wire  data  have 
been  obtained  at  eleven  chordwise  stations  along  the  airfoil,  from  x/c=0.13  to  x/c=0.98.  All  data 
have  been  obtained  at  a  fixed  z  location  of  0.5"  above  the  airfoil  shoulder  (z/c=0.04).  At  each 
location  data  have  been  obtained  at  a  sampling  rate  of  10  kHz,  in  Ik  4k  or  16  k  blocks.  The 
present  results  are  for  1  k  data  records,  corresponding  to  0.1  sec  time  records.  This  has  been 
done  to  minimize  the  time  required  to  obtain  a  PDV  data  set,  in  order  to  minimize  time  for  iodine 
cell  stem  temperatures  to  drift.  Such  drift  of  the  temperatures  of  the  iodine  cell  stems  has  been 
identified  as  the  key  source  of  the  apparently  randomly  varying  offset  errors  in  mean  velocity 
results.  From  these  1  k  time  series,  mean  and  RMS  velocities,  as  well  as  time  autocorrelation 
coefficients  and  power  spectra  have  been  computed  for  both  the  PDV  and  hot  wire  data,  to  assess 
the  accuracy  of  the  PDV  system  for  turbulence  measurements.  Also,  similar  data  are  presently 
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being  acquired  in  a  series  of  one  inch  diameter  jets:  a  uniform  jet,  an  annular  jet,  and  a  swirling 
jet. 

7.2  Results 

Early  PDV  data  repeatability,  as  documented  in  the  thesis  by  Ramanath  (1997),  was  poor. 
However,  improved  cell  calibration  procedures  (James,  1997)  have  significantly  increased  PDV 
and  DGV  system  accuracy.  Both  single  and  2-component  PDV  data  measured  on  a  rotating 
wheel  have  been  presented  by  James  (1997),  and  by  Kuhlman,  et  al.  (1997).  Total  wheel 
velocity  range  for  these  measurements  was  57  m/sec,  so  the  observed  velocity  error  magnitudes 
of  approximately  ±0.6- 1.2  m/sec,  correspond  to  1-2  %  errors,  which  is  quite  good.  Also,  the 
standard  deviations  of  the  actual  PDV  wheel  velocity  data  points  from  the  least  squares  linear 
curve  fits  have  been  found  to  be  0.5-0.7  %. 

Two-component  PDV  data  obtained  from  a  traverse  across  the  exit  of  the  fully-developed 
pipe  flow  at  a  nominal  Reynolds  number  of  76,000,  have  been  previously  presented  by  Kuhlman 
(1998).  The  axial  mean  velocities  agree  well  with  results  from  a  pitot-static  probe  survey,  but 
circumferential  mean  velocities  display  an  offset  error  on  the  order  of  2-5  m/sec.  Turbulent 
velocity  levels  agree  with  hot  wire  data  of  Laufer  (1954).  Significant  difficulties  occur  near  the 
pipe  walls,  both  due  to  reduced  signal-to-noise  levels  due  to  less  smoke,  and  to  secondary 
scattering  and  reflections  of  the  scattered  light  off  of  the  pipe  walls.  However,  the  greatest 
difficulty  was  in  obtaining  accurate,  reliable  zero  velocity  values. 

Recent  PDV  data  has  been  compared  with  hot  wire  data  by  Kuhlman  and  Webb  (1999) 
for  flow  over  an  NACA  0012  airfoil  model  at  zero  degrees  angle  of  attack.  PDV  streamwise 
mean  velocities  show  a  total  variability  of  about  5  to  6  m/sec  for  the  4  runs;  this  variability  is 
consistent  with  earlier  PDV  data  (Kuhlman,  1998)  and  DGV  data  (Naylor  and  Kuhlman,  1999). 
PDV  RMS  velocities  agree  with  corresponding  hot  wire  data  to  within  about  0.5m/sec,  and 
display  a  similar  level  of  repeatability  from  run  to  run.  The  computed  PDV  correlation 
coefficients  agree  well  with  hot  wire  correlation  coefficients;  both  autocorrelations  decay  in  less 
than  10  msec.  However,  PDV  power  spectra  show  a  larger  amount  of  noise  at  frequencies  above 
about  1000  Hz. 
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Recent  2-component  PDV  data  in  a  circular  jet  are  shown  in  Figs.  7.1  and  7.2.  Mean  and 
RMS  data  at  the  jet  exit  (x/D=0.25,  where  D=exit  diameter)  are  shown  in  Fig.  7.1,  while  similar 
results  are  shown  in  Fig.  7.2  for  x/D=6.  Mean  streamwise  velocities  at  the  exit  (Fig.  7.1)  again 
show  about  5  m/sec  of  offset  from  preliminary  pitot  probe  results  on  the  centerline  at  the  exit. 
Measured  PDV  mean  circumferential  velocity  varies  by  about  2  m/sec  from  the  correct  value  of 
zero  m/sec.  RMS  velocities  are  between  1  and  2  m/sec  on  the  centerline  (Fig.  7.1).  At  x/D=6, 
the  PDV  mean  streamwise  velocity  profile  (Fig.  7.2)  appears  Gaussian  in  shape,  as  expected. 
Both  the  mean  streamwise  and  circumferential  velocities  display  offset  errors  of  4-5  m/sec  from 
the  correct  values.  RMS  velocities  have  risen  to  4-10  m/sec.  Detailed  comparisons  are  now 
being  made  between  PDV  and  hot  wire  data  for  this  turbulent  circular  jet  flow  field.  Also,  the 
same  jet  flow  facility  has  been  modified  to  allow  PDV  and  hot  wire  data  to  be  obtained  for  a 
swirling  jet  as  well  as  for  an  annular  jet. 

Recent  analysis  of  the  offset  error  has  indicated  that  it  is  largely  due  to  the  random, 
uncorrelated  variations  in  iodine  cell  stem  temperatures.  This  stem  temperature  variation  has 
been  observed  to  vary  with  a  short  term  RMS  of  0.1  degrees  C  (Naylor,  1998).  This  has  been 
found  to  correspond  to  an  error  in  the  computed  Doppler  shift  frequency  of  7  MHz,  which  then 
corresponds  to  a  mean  velocity  error  of  from  2  to  10  m/sec,  depending  on  the  geometry  and 
viewing  direction  of  the  PDV  system.  This  level  of  velocity  error  is  consistent  with  the  observed 
mean  velocity  offset  error. 

A  recent  upgrade  has  been  performed  of  the  iodine  cells  and  beamsplitters  used  in  the 
present  system.  Vapor-limited  iodine  cells  have  recently  been  installed,  along  with  polarization- 
insensitive  beamsplitters.  Now,  both  PDV  and  DGV  data  are  being  taken  in  the  one-inch 
diameter  circular  jet  facility,  for  comparison  with  pitot-static  and  hot  wire  data  that  has  recently 
been  obtained.  It  is  expected  that  the  new,  vapor-limited  cells  will  largely  eliminate  the  mean 
velocity  offset  error  that  has  been  observed  in  all  point  and  DGV  data  taken  to  date.  Also,  a 
DGV  system  using  line  scan  cameras  is  under  development,  and  a  higher  power,  pulsed  YAG 
laser  will  be  added  to  the  existing  system.  All  of  these  efforts  are  aimed  at  both  improving 
measurement  accuracy  and  allowing  turbulence  measurements  to  be  made  with  the  improved 
DGV  and  PDV  systems. 
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Chapter  8  Turbulence  Predictions  in  1C  Engines 


In  this  section  we  present  our  most  important  findings  from  application  of  the  LES 
technique  to  typical  engine  cylinder  geometries.  First,  the  flow  dynamics  due  to  piston  motion 
during  the  compression  and  expansion  strokes  is  studied.  Here  the  primary  instability  is  induced 
by  the  presence  of  a  bowl,  and  the  squish  flow  combined  with  an  imposed  swirl  component.  The 
initial  conditions  for  these  simulations  were  only  approximate  and  are  not  reminiscent  to  the 
intake  stroke.  Next,  the  flow  inside  the  cylinder  during  the  intake  stroke  and  compression  strokes 
without  combustion  was  simulated  including  the  valve  dynamics.  Here,  the  objective  was  to 
investigate  the  generation  of  turbulence  during  the  intake  and  subsequently  its  decay  during 
compression  in  addition  to  demonstrating  the  predictive  capability  of  the  LES  technique  when 
applied  to  in-cylinder  turbulence.  Finally,  two  cases  were  simulated  with  combustion  (1)  was  the 
continuation  of  the  intake  flow  case  with  spray  injection  and  subsequent  combustion,  (2) 
combustion  with  a  piston  bowl,  again  without  the  influence  of  the  turbulence  generated  during 
the  intake  stroke.  The  results  are  compared  to  available  experiments  whenever  possible. 

8.1  Bowl  Induced  Flow  Instability 

This  study  (Yavuz  et  al.,  1998)  focuses  on  capturing  flow  instability,  and  possibly 
turbulence,  induced  by  a  typical  piston  bowl  geometry  by  simulating  the  flow  dynamics  using 
the  KIVA-3  code.  Both  two-  and  three-dimensional  simulations  were  performed  using  relatively 
fine  grid  resolution  in  a  piston-bowl  assembly  with  various  wedge  angles  in  the  circumferential 
direction.  The  differences  in  solutions  using  a  subgrid  scale  (SGS)  model  are  elucidated. 
Simulations  with  no  turbulence  model  seem  to  capture  a  high  frequency  of  the  instabilities  with 
significant  amplitudes.  The  r.m.s.  (root  mean  square)  values  of  the  fluctuations  have  been 
calculated  at  selected  points  inside  the  cylinder.  An  assessment  of  the  influence  of  numerical 
parameters  on  the  predictions  is  made. 

The  results  are  not  compared  with  any  experimental  work  directly  because  such  a 
comparison  requires  strict  matching  of  the  intake  flow  situation,  which  is  not  done  in  the  present 
study.  Moreover,  the  experiments  most  of  the  time  present  ensemble  averages  which  include 
cycle-to-cycle  variations,  and  they  can  not  be  easily  isolated.  However,  the  experiments  by 


8S 


Catania  et  al.,  (1995),  and  Valentino  et  al.,  (1998)  are  used  as  a  guide  to  make  qualitative 
assessment  of  the  predictions.  The  simulations  are  enlightening  in  that  the  turbulence  inside  an 
engine  cylinder  can  be  predicted,  and  further  the  influence  of  cylinder  bowl  geometry  can  be 
isolated  from  that  of  the  turbulence  generated  by  other  mechanisms  such  as  shear  layers  and  the 
swirl  induced  during  the  intake  stroke. 

8.1.1  Application 

The  flow  field  inside  a  typical  piston-cylinder  assembly  (Figure  8.1,  2-D,  axisymmetric) 
was  simulated  under  motored  conditions  (i.e.,  without  combustion)  at  two  speeds,  namely  600 
RPM  and  1500  RPM.  Only  the  compression  and  expansion  strokes  are  simulated;  cycle-to-cycle 
variations  are  not  considered.  These  simulations  were  performed  using  a  maximum  time  step  of 
0.15  ps  over  a  compression  and  expansion  stroke  repeated  twice.  The  focus  is  on  the  prediction 
of  flow  instabilities  induced  by  the  bowl  geometry  and  the  squish  flow,  which  may  contribute 
significantly  to  the  overall  turbulence  generation.  The  numerical  mesh  used  in  the  simulations  is 
depicted  in  Figure  8.2.  Two  numerical  meshes  of  20,000  and  130,000  vertices  were  used  with  a 
wedge  angle  of  0.5  and  180  degrees  in  the  circumferential  direction,  respectively.  The  maximum 
cell  size  for  the  wedge  calculations  was  0.5x0.9  mm  in  the  radial  and  axial  direction.  The  3-D 
simulation  with  180-degree  wedge  angle  span  was  investigated  with  18  cells  in  the 
cirumferential  direction  to  see  the  influence  of  the  domain  size.  For  this  case  the  median  cell  size 
in  the  tangential  direction  was  3.0  mm.  As  reported  by  Celik  and  Yavuz  (1997),  the  minimum 
integral  length  scale  occurs  at  the  end  of  the  compression  stroke  and  is  in  the  order  of  a  few 
millimeters  for  automotive  size  engines. 

Early  calculations  have  revealed  that  it  was  possible  to  capture  flow  instability  with  the 
present  grid  resolution  even  at  the  lower  engine  speed  of  600  RPM.  However,  the  selection  of  the 
time  step  played  a  crucial  role  in  the  resulting  amplitude  and  the  frequency  of  oscillations.  It  was 
possible  to  totally  destroy  the  flow  instability  (most  probably  induced  by  the  flow  separation  at 
the  comer  where  the  bowl  meets  the  piston  head)  either  by  increasing  the  time  step,  or  by  simply 
using  the  standard  k-e  turbulence  model,  which  is  known  to  be  overly  diffusive.  The  time  steps 
used  in  the  results  presented  in  this  paper  were  selected  carefully  after  much  experimentation. 
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Figure  8.1  2-D  Streamlines  at  CA  210°  ATDC  (with  mirror  image) 


Figure  8.2  Computational  mesh  for  2-D  simulation  and  specific  engine  data 

As  it  was  mentioned  above  the  quasi-second  order  upwind  (QSOU)  scheme  was  used  for 
the  convective  terms,  and  central  differencing  for  the  diffusion  terms.  Normally  the  central 
differencing  (CD)  scheme  should  be  used  for  the  momentum  equations.  However,  CD  is  usually 
an  unstable  scheme  when  used  for  convection.  The  stability  can  be  secured  by  adding  some 
diffusion,  via  either  SGS  model  or  artificially.  Therefore,  prediction  of  turbulent  fluctuations  via 
the  QSOU  scheme,  which  is  strictly  monotonic,  was  preferred  at  the  expense  of  some  loss  of 
formal  accuracy  near  sharp  gradients  that  may  be  present  in  the  flow  field. 

The  simulations  were  started  at  CA  90°  BTDC.  As  an  initial  condition  the  radial  velocity 
was  taken  as  zero  and  the  axial  velocity  varied  linearly  from  piston  velocity  at  the  piston  face  to 
zero  at  the  cylinder  head.  The  tangential  velocity  was  defined  by  the  Bessel  function  dependent 


on  engine  speed  as  provided  in  the  code.  Moreover,  periodic  boundary  conditions  were  applied 
on  the  front  and  back  face  of  the  wedge,  whereas  the  law  of  the  wall  was  applied  on  the  wails. 
And  in  the  case  the  SGS  model  was  applied,  an  initial  turbulence  intensity  of  10%  and  Lsos  =  0.4 
cm  was  used,  whereas  without  the  SGS  model  no  turbulence  quantities  were  initiated  at  all. 

Cycle  resolved  phase  averaging  was  applied  to  the  calculated  data  to  separate  the  mean 
and  fluctuating  components  of  the  instantaneous  velocity  as  explained  in  Yavuz  et  al.  (1998). 

8.1.2  Results 

An  example  of  the  complexity  of  the  flow  field  as  predicted  with  the  present  model  is 
shown  in  Figure  8.1.  This  figure  is  indicative  of  the  level  of  details,  which  can  be  captured  by  the 
present  simulations.  The  variation  of  the  instantaneous  velocity  components  with  time  at  location 
A  (see  Fig.  8.2),  which  is  fixed  with  respect  to  the  bowl,  1.0  mm  below  the  piston  face  and  10.6 
mm  from  the  axis,  are  depicted  along  with  the  time  averaged  values  in  Figure  8.3  for  the 
compression-expansion  stroke.  Here,  the  two  indicated  TDC's  (Figure  8.3  and  8.4)  are  repetitions 
of  the  compression/expansion  stroke. 

As  it  was  expected,  the  fluctuations  in  the  tangential  direction  are  much  smaller 
compared  to  the  other  two  directions.  This  is  a  direct  consequence  of  the  pseudo  two- 
dimensional  nature  of  the  present  calculations  with  0.5  degree  wedge  angle.  Both  the  large-scale 
time  variations  of  the  velocity  components  induced  by  the  piston  motion,  as  well  as  the  random 
"turbulent"  fluctuations  are  captured.  Figure  8.4  shows  these  random  velocity  fluctuations  as  a 
function  of  the  crank  angle.  It  is  seen  that  the  turbulence  intensity  is  somewhat  reduced  during 
the  second  compression/expansion  stroke;  however,  it  persists  at  a  considerable  level  indicating 
that  the  numerical  diffusion  is  sufficiently  small  not  to  render  turbulent  dynamics  of  the  flow. 
Whether  these  fluctuations  will  continue  to  persist  over  several  cycles  remains  to  be  seen.  Figure 
8.4  and  8.5  indicate  that  there  are  significant  changes  occurring  near  the  TDC  both  in  the 
magnitude  and  frequency  of  the  velocity  fluctuations.  As  it  is  shown  in  Figures  8.5  and  8.6  the 
time  averaging  interval,  as  well  as  the  time  step  used  for  the  numerical  simulation  has  significant 
influence  on  the  intensity  of  the  fluctuations. 
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Figure  8.3  Instantaneous  velocity  components  at  point  A  (without  SGS  model) :  pseudo  2-D 
calculations  for  a  motored  engine  at  1 500  RPM 


Figure  8.4  Instantaneous  fluctuating  velocity  components  at  point  A  (without  SGS  model): 
pseudo  2-D  calculations  for  a  motored  engine  at  1500  RPM 


Figure  8.5  Influence  of  time  averaging  interval  on  various  flow  quantities:  pseudo  2-D 
calculations  for  a  motored  engine  at  1500  RPM 
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It  should  be  pointed  out  that  the  predicted  "turbulent"  fluctuations  are  significantly 
affected  by  the  time  step  (Figure  8.6a).  The  smaller  the  time  steps,  the  better  is  the  resolution  of 
turbulence  scale,  and  the  higher  is  the  turbulent  intensity.  Figure  8.6b  shows  the  influence  of 
SGS  turbulence  model,  which  essentially  has  a  damping  effect.  At  this  grid  resolution,  which  is 
relatively  coarse  for  LES  simulations  we  believe  that  the  results  without  the  SGS  model  should 
be  better  than  the  ones  with  a  SGS  model. 


(c) 
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Figure  8.6  Influence  of  various  parameters  on  calculated  quantities:  pseudo  2-D  calculations 
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The  predicted  fluctuations  are  anisotropic  (Fig.  8.6)  even  very  near  the  TDC.  The 
magnitudes  of  the  fluctuations  are  in  the  range  0.5- 1.0  m/sec  in  radial,  1. 0-2.0  m/sec  in  the  axial 
and  0.5 -1.5  m/sec  in  the  tangential  directions,  respectively.  For  a  similar  engine  configuration, 
Catania  and  Spessa  (1996)  measured  velocity  fluctuations  in  the  order  of  10  m/s  for  the  cycle- 
resolved  turbulence  intensities,  and  about  5m/s  for  the  r.m.s.  fluctuations  of  the  in-cyclinder 
mean  velocity.  According  to  Catania  (1998),  the  present  computations  should  be  compared  to 
their  cycle-resolved  turbulence  intensities.  In  another  experimental  study  Auriemma  et  al., 
(1998)  reported  measured  rms  velocity  fluctuations  during  the  compression/expansion  stroke  of  a 
motored  engine  at  1500  RPM.  These  indicated  maximum  values  of  2  m/s  for  the  radial  and  3.5 
m/s  for  the  tangential  directions  respectively.  The  measured  values  included  both  cycle-to-cycle 
fluctuations  and  in-cycle  fluctuations.  Valentino  et  al.,  (1998)  presented  a  method  to  reduce  the 
effect  of  cycle-to-cycle  variations  of  mean  motion  on  turbulence.  Their  results  seem  to  indicate 
that  the  cycle-resolved  turbulence  can  be  augmented  by  as  much  as  30%  due  to  cycle-to-cycle- 
variations.  In  any  case,  the  experimental  data  indicate  that  the  predicted  turbulence  intensities  are 
much  lower  compared  to  those  measured,  although  direct  comparison  is  not  possible  at  this 
stage.  This  may  be  mainly  due  to  the  fact  that  reminiscent  turbulence  generated  during  the  intake 
stroke  was  not  accounted  for. 


Figure  8.7  Streamlines  to  show  the  influence  of  three-dimesionality  on  the  flow  field  at  CA  18° 

ATDC  at  a)  0°  and  b)  90°  wedge  angles 


In  Figure  8.6c,  the  influence  of  engine  speed  is  depicted.  There  is  a  significant  increase  in 
predicted  turbulent  intensity  with  increasing  engine  speed  as  has  been  observed  in  experiments 
(see  e.g.,  Corcione  and  Valentino,  1994;  Catania  and  Spessa,  1996). 

The  3-D,  180°  wedge  calculations  with  130000  vertices  did  not  show  any  significant 
change  in  the  predicted  rms  fluctuations.  This  could  be  due  to  the  stabilizing  effect  of  the 
arbitrarily  imposed  swirl  velocity  profile  in  addition  to  the  coarse  grid  used  in  the  tangential 
direction.  However,  the  streamlines  shown  on  two  perpendicular  planes  of  the  3-D  geometry  (in 
Figures  8.7a  and  b),  reveals  some  expected  3-D  effects. 

8.1.3  Discussions 

Three  dimensional,  and  pseudo-two  dimensional  (i.e.,  a  small  wedge  angle),  time 
accurate  calculations  have  been  performed  to  predict  the  flow  field  inside  a  typical  diesel  engine 
piston-bowl  assembly.  The  engine  was  run  under  motored  conditions  at  600  and  1500  RPM, 
spanning  only  the  compression  and  expansion  strokes.  The  objective  was  to  capture  the  flow 
instabilities,  and  possibly  turbulence  induced  by  the  bowl  geometry  alone,  i.e.  without  the 
influence  of  residual  turbulence  generated  during  the  intake  stroke.  It  was  shown  that  this  was 
possible  by  using  the  KIVA-3  code  with  a  fine  grid  resolution  and  relatively  small  time  steps. 
The  bowl  induced  “turbulence”  seems  to  be  persistent  over  a  few  cycles  without  recharging 
induced  by  the  intake  stroke.  The  small  scale  flow  instabilities  and  unsteadiness  observed  in  the 
computer  simulations  seem  to  be  the  results  of  imsteady  flow  separation  that  occurs  at  the  sharp 
comer  of  the  piston  bowl  and  the  flat  piston  head.  This  is  in  accordance  with  the  combustion 
improvements  that  are  observed  experimentally  when  re-entrant  type  bowl  geometry  is  used.  The 
magnitude  of  the  calculated  r.m.s.  velocity  fluctuations  seem  to  be  much  lower  than  measured 
ones  under  similar  geometry  and  operating  conditions.  But  the  calculations  do  not  include  the 
significant  contribution  that  arises  from  the  turbulence  generated  during  the  intake  stroke,  nor  do 
they  include  the  significant  cycle-to-cycle  variations.  Although  the  grid  resolution  used  in  the 
present  study  may  be  considered  as  relatively  coarse  for  LES,  this  study  demonstrated  that  there 
is  a  potential  to  predict  in-cylinder  turbulence  for  IC  engines  using  large  eddy  simulations  with 
relatively  modest  computational  resources  at  the  workstation  level. 
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8.2  Prediction  of  In-Cylinder  Turbulence  for  1C  Engines 


Results  of  under-resolved  large  eddy  simulations  of  an  intake  and  a  compression  stroke  in 
a  typical  IC  engine  geometry  are  presented.  A  modified  KIVA-3V  program  is  used  with  a  special 
precaution  to  keep  the  time  accuracy  at  a  level  acceptable  for  LES.  A  Smagorinsky  sub-grid 
scale  model  is  used  with  improved  spatial  discretization  scheme.  The  simulations  could  achieve 
a  spatial  resolution  of  turbulent  eddies  on  the  order  of  several  millimeters  and  temporal 
resolution  of  velocity  fluctuations  on  the  scale  of  lO'^s.  The  predicted  growth  and  the  subsequent 
decay  of  turbulence  during  the  intake  phase  agree  well  with  experiments.  Statistics  of  turbulent 
velocity  fluctuations  is  analyzed  and  compared  with  experimental  studies.  Combustion  of  n- 
Tetradecane  is  modeled  using  a  7-reaction  Arrhenius  mechanism  and  Magnussen  eddy- 
dissipation  combustion  model.  A  standard  KIVA  spray  model  is  used  to  model  injection  and 
evaporation  of  the  fuel.  Auto-ignition  of  the  gaseous  phase  of  the  fuel  occurs  close  to  10“  CA 
BTDC.  The  influence  of  combustion  on  turbulence  is  analyzed. 

8.2.1  Numerics 

Modifications  were  made  to  the  KIVA  code  to  allow  a  fairly  simple  Smagorinsky  SGS 
model  where  the  eddy  viscosity  is  calculated  from  the  algebraic  relation  (Eqn.  1 .6).  The  model 
constant  Cs  was  set  equal  to  0.2,  a  typical  value  used  in  the  literature  (see  e.g.  Rodi  et  al.,  1997). 
The  sub-grid  length-scale  was  estimated  as  an  average  computational  cell  dimension  (see  Smith 
etal.,  1998). 

In  some  of  the  cases  presented  “no  SGS  turbulence  model”  was  used.  This  requires,  in 
most  instances,  some  numerical  diffusion  for  stability,  which  is  accomplished  via  the  use  of 
QSOU  (Quasi  Second  Order  Upwind)  scheme  instead  of  the  central  differencing  (CD)  which  has 
no  diffusion  error,  or  by  using  a  combination  of  CD  and  upwind  differencing  heavily  biased 
towards  CD. 

The  above  turbulence  models  were  employed  in  conjunction  with  the  commonly  used 
law-of-the-wall  boundary  condition  which  is  implemented  in  KIVA-3.  The  basic  assumption 
here  is  that  the  interaction  between  the  modeled  near  wall  region  and  the  outer  region  is  weak  as 
observed  experimentally  by  Brooke  and  Hanratty,  (1993). 
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Convective  terms  are  advanced  explicitly  in  time  and  the  diffusion  terms  are  advanced 
explicitly,  implicitly,  or  semi-implicitly.  The  degree  to  which  the  diffusion  terms  are  implicitly 
discretized  is  based  on  a  combination  of  stability  and  efficiency  considerations.  Sub-time  steps 
(referred  to  as  sub-cycles)  are  taken  to  advance  the  convective  terms.  The  time  step  in  each  sub¬ 
cycle  is  based  on  Courant  stability  considerations.  The  convective  terms  are  advanced  with  time 
steps  that  are  the  same  or  smaller  than  what  is  used  for  the  diffusion  terms.  Though  the  overall 
time  accuracy  for  convection  terms  is  only  of  the  first-order,  the  global  time  step  is  based  on 
several  considerations  including  stability  and  several  accuracy  constraints.  To  increase  the  time 
accuracy  even  further  the  time  step  was  confined  to  the  interval  [1.5T0'^  -  5T0'’  sec]  far  below 
the  Kolmogorov  time  scale,  which  is  in  the  order  of  10*^  -10'^  seconds  for  a  typical  automotive 
size  engine  (Celik  and  Yavuz,  1997).  This  precaution  compensates  for  the  first  order  time 
accuracy  of  convective  terms  and  guarantees  the  time  resolution  required  for  LES. 

8.2.2  Engine  Application 

The  computational  procedure  described  above  was  applied  to  a  typical  engine  geometry. 
The  cylinder  assembly  was  prototyped  after  the  Diesel  engine  used  in  the  experiments  of  Catania 
et  al.  (1996)  with  the  following  dimensions:  bore  -  79.5mm,  stroke  -  86mm,  compression  ratio  - 
18,  inlet  valve  opening  -  5  deg  BTDC,  inlet  valve  closure  -  55  deg  ATDC,  exhaust  valve  opening 
-  55  deg  BTDC,  exhaust  valve  closure  5  deg  ATDC.  Piston  speed  was  set  to  1600  rpm.  A 
realistic  valve  geometry  with  the  swirling  intake  configuration  was  used  (Figure  8.8).  The 
computations  were  performed  on  two  grids  consisting  of  220,000  and  440,000  nodes  including 
the  intake/exhaust  ducts.  The  average  grid  size  was  in  the  order  of  one  millimeter.  A  curvilinear 
block-structured  coordinate  system  was  used  in  preference  to  the  cylindrical  coordinate  system 
because  the  former  avoids  high  grid  distortions  at  the  cylinder  axis.  The  first  set  of  computations 
on  both  the  coarse  and  the  fine  grids  was  carried  out  without  any  turbulence  model  and  with  the 
QSOU  numerical  scheme.  Later  the  computations  on  the  fine  grid  were  repeated  with  the 
Smagorinsky  sub-grid-scale  turbulence  model  and  the  second  order  discretization  scheme  for  the 
convective  terms  (as  described  by  Smith  et  al.,  1999).  The  exhaust  duct  pressure  was  set  equal  to 
the  atmospheric  pressure  and  the  intake  pressure  was  set  40%  above  the  former,  to  prevent  the 
in-cylinder  gas  from  entering  the  intake  duct  during  the  first  phase  of  intake.  No  random  forcing 
was  used  at  the  inlet.  It  should  be  noted  that  the  simulated  engine  configuration  and  the  one 
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studied  experimentally  by  Catania  et  al.  are  not  exactly  the  same.  The  simulated  engine  has  a  flat 
piston  head,  whereas  the  experimental  one  had  a  specifically  designed  piston  bowl,  which  is  off- 
centered  with  respect  to  the  cylinder  axis.  These  differences  were  due  to  the  limitations  of  the 
computer  code  and  the  currently  available  computational  resources.  The  intake  valve  geometry 
and  movement  were  modeled  in  accordance  with  the  real  engine  parameters,  but  a  generic  profile 
was  used  to  describe  the  shape  of  the  upper  portions  of  the  valves  (Figure  8.8b). 


Figure  8.8  (a)  Geometry  of  the  simulated  engine;  (b)  Velocity  distribution  in  the  valve  region 


The  Smagorinsky  sub-grid-scale  model  and  the  Magnussen  Eddy-Dissipation  model 
(Magnussen,  1976)  were  used  to  account  for  the  effects  of  turbulence  and  combustion 
respectively.  Time  step  was  varied  depending  on  the  time-scale  of  the  limiting  process.  Thus,  in 
the  evaporation  phase  (285-360  CA)  the  time  scale  could  reach  the  order  of  lO’^s,  which  was  due 
to  the  limits  imposed  by  the  evaporation  process  of  smallest  droplets.  A  standard  KIVA  spray 
model  is  used  in  all  the  simulations  to  account  for  injection,  initial  droplet  size  distribution  and 
evaporation  of  the  fuel.  The  initial  gas  and  wall  temperatures  of  450°K  were  used.  The  average 
droplet  diameter  of  the  spray  was  set  equal  to  lO'^cm  and  STO'^cm  for  the  case  with  a  bowl  and 
two-valve  assembly  respectively. 


8.2.3  Statistical  Data  Processing 

The  fluctuating  component  u(t)  was  obtained  by  subtracting  the  mean  value,  U(t) ,  from 
the  instantaneous  velocity 
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u(t)  =  U(t)-U(t) 


(8.1) 


The  mean  velocity  was  computed  by  continuous  time  averaging  (phase  averaging)  over 
the  interval  T  corresponding  to  the  crank  angle  interval  d)  by  discrete  integration  corresponding 
to. 

_  I  T/2  ,  0/2 

U(t)  =  -  |U(t  +  T)dT  =  —  JU(0(t)  +  &)d&  (8.2) 

T_t/2  *^-0/2 

The  auto  correlation  function  R*uu  and  the  power  spectral  density  were  computed 
according  to  the  following  formulae 

Ruu(^)  =  ^  ju  (t-T/2)u(t  +  T/2)dt  (8.3) 

a  ^70-1/2 
T/2 

Euu(0  =  4  I  Ruu(T)w(T)cos(27ifT)dT  (8.4) 

0 

where  To  corresponds  to  the  crank  angle  of  Oq,  which  is  the  middle  of  the  integration 
interval.  The  Hanning  window  function  w(t)  was  applied  to  decrease  the  variance  in  the 
spectrum.  A  similar  technique  was  applied  in  the  experimental  study  of  Catania  et  al  (1996), 
which  was  used  for  comparison.  The  difference  between  our  definitions  of  E*uu.  R*uu  and  those 
used  by  Catania  et  al.,  Euu,  Ruu  is  that  in  the  latter  case  the  ensemble  averaging  over  cycles  is  also 
applied  in  addition  to  phase  averaging. 

8.2.4  Results  of  Intake  Turbulence  Predictions: 

The  contours  of  the  magnitude  of  the  velocity  vector  depicted  in  Figure  8.9  illustrate  the 
details  of  the  three-dimensional  flow  structures  at  two  different  crank  angles  that  are  captured  by 
the  simulations.  The  figure  shows  a  largely  asymmetric  armular  intake  jet  with  significant 
number  of  eddies  forming  in  the  shear  layers  surrounding  the  main  body  of  the  jet.  The  sizes  of 
the  smallest  eddies  captured  are  proportional  to  the  grid  size  (~1  mm).  The  integral  length  scale 
is  estimated  (Celik  and  Yavuz,  1997)  to  attain  a  minimum  value  of  a  few  millimeters  at  the  TDC 
proportional  to  the  squish/clearance  height.  The  details  of  the  flow  inside  the  intake  duct  were 
smeared  (filtered)  out  by  the  relatively  coarse  grid  in  that  region. 
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Figure  8.10  shows  the  fluctuating  velocity  components  and  its  fluctuating  part  located  at 
a  point  approximately  5  mm  below  the  cylinder  head  in  the  middle  between  the  centerlines  of  the 
intake  and  exhaust  ducts.  As  such  this  particular  point  can  be  compared  with  one  measurement 
point  of  Catania  and  Spessa  (1996)  which  is  5  mm  below  the  cylinder  head,  and  close  to  the 
cylinder  axis.  The  reader  is  reminded  that  the  geometry  of  Catania's  experiment  is  different  from 
the  computational  model  in  that  there  is  no  cylinder  bowl  in  the  latter  case.  Because  of  that  and  a 
slight  difference  in  special  locations  where  the  measured  and  computed  velocities  were  retrieved 
it  was  not  possible  to  predict  the  mean  flow-field  exactly.  But  our  main  concern  was  to 
reproduce  the  trend  in  the  fluctuating  component,  which  is  seen  to  agree  well  with  experimental 
data  (Figure  8.10b).  The  predicted  cycle  resolved  mean  velocity  (not  depicted  here)  showed 
roughly  the  same  trend  as  in  experiments  but  with  a  slightly  lower  magnitude.  It  is  noteworthy  to 
see  that  the  velocity  fluctuations  decrease  significantly  towards  the  BDC  (-180  CA  degree)  as 
observed  in  the  measurements. 


105°  crank  angle. 


120°  crank  angle. 

Figure  8.9  Absolute  value  of  the  velocity  vectors  (440000  node  case). 
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One  of  the  early  experimental  study  demonstrating  the  decay  of  turbulence  during  the 
expansion  stroke  is  that  of  Semenov  and  Sokolik,  (1958).  In  these  experiments  the  piston  head 
was  flat  which  corresponds  to  the  computational  setup,  but  the  intake  valve  was  located  co¬ 
axially  relative  to  the  cylinder.  The  hot  wire  anemometer  used  was  insensitive  to  the  direction  of 
velocity  thus  only  the  instantaneous  magnitude  of  velocity  was  measured.  According  to  the 
authors  their  measurements  underestimate  the  amplitude  of  velocity  fluctuations  by  15-20%.  In 
Figure  8.1 1  the  measured  and  computed  time  dependence  of  the  absolute  value  of  velocity  in  the 
plane  normal  to  the  cylinder  axis  is  shown.  The  general  trend  of  decreasing  velocity  variations 
(both  in  the  mean  and  fluctuating  components)  can  be  observed  in  both  cases.  As  our  grid 
dependence  study  showed  doubling  the  number  of  nodes  in  the  grid  revealed  more  fine  structures 
in  the  flow  and  led  to  higher  amplitudes  of  velocity  fluctuations.  Nevertheless,  the  general  trend 
of  turbulence  decay  was  unchanged  in  both  cases  (Figure  8.12).  It  is  important  to  note  that  these 
particular  measurements  (both  Catania  et  al.’s  and  Semenov’s  data)  show  cycle-resolved 
turbulence  rather  than  ensemble  averages. 

Finally  we  present  the  normalized  autospectral  density  function  (i.e.  energy  spectra)  in 
comparison  with  experiments  in  Figure  8.13  at  a  representative  point  which  roughly  corresponds 
to  the  measurement  point.  The  agreement  of  both  the  energy  content  and  the  frequency  ranges 
captured  in  the  present  simulations  are  good  in  spite  of  the  differences  in  the  cylinder 
configurations  as  explained  above.  The  calculated  spectra  looks  wavy  as  compared  to  the 
experiments,  but  this  is  most  probably  due  to  the  smoothing  effect  of  the  double  averaging 
(including  ensemble  averages)  applied  to  measure  data  (see  Equation  11,  Catania  et  al  (1996)). 
Indeed  when  the  average  of  the  wavy  autocorrelation  function  is  used  to  compute  the  energy 
spectra  the  wiggles  seen  in  the  computed  curve  reduced  drastically.  Figure  8.13  is  seen  as  a 
confirmation  of  the  fact  that  the  present  methodology  is  capable  of  resolving  a  significant  portion 
of  the  inertial  sub-range  in  the  energy  spectra. 
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Figure  8. 1 3  Power  density  spectra:  a)  computed  (engine  without  a  bowl),  b)  measured  (engine 
with  a  bowl):  r  =  radial,  s  =  tangential,  z  =  axial  wire 

8.2.5  Combustion  Results: 

Two  types  of  simulations  were  performed  to  approach  a  realistic  case  of  three- 
dimensional  fine-grid  combustion  simulation.  At  the  first  stage  a  2D  grid  with  71x75  nodes  in 
the  cylinder  region  and  31x23  nodes  in  the  bowl  region  was  used.  In  3D  simulations  this  grid 
was  extended  to  150,000  nodes  by  a  simple  revolution  of  the  2D  grid.  A  7-reaction  reduced 
mechanism  for  n-Tetradecane'  with  Arrhenius  rate  expressions  was  used  at  this  stage.  A  standard 
KIVA  spray  model  is  used  in  all  the  simulations  to  account  for  injection,  initial  droplet  size 
distribution  and  evaporation  of  the  fuel.  The  initial  gas  and  wall  temperatures  of  450°  K  were 
used.  The  average  droplet  diameter  of  the  spray  was  set  equal  to  ITO"^  cm  and  3-10'^  cm  for  the 
case  with  a  bowl  and  two-valve  assembly,  respectively.  No  SGS  model  was  used  for  the  auto¬ 
ignition  of  the  fuel  in  its  gaseous  phase,  which  occurs  near  TDC  approximately  at  352°  CA. 
Figure  8.14  shows  the  temperature  contours  of  this  case  on  a  vertical  cross-section  just  after 
TDC.  The  perfectly  symmetric  initial  conditions  were  broken  during  the  three-dimensional 
transient  calculations,  and  resulted  in  a  three  dimensional  skewed  flame  structure.  In  the  later 
case  with  flat  piston,  the  asymmetry  was  already  in  the  flow  due  to  valve  arrangements.  Here  the 
combustion  mostly  occurs  near  the  wall.  The  predicted  flame  front  is  rather  smooth.  Due  to  the 
very  simple  combustion  model  used,  the  wrinkled  laminar  flamelet  structure  could  not  be 
captured.  Nevertheless  these  calculations  showed  that  it  is  feasible  to  do  engineering  LES 
combustion  in  IC  engines  using  the  widely  used  engine  code  KIVA.  Work  is  underway  to  refine 
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the  combustion  sub-grid  scale  models  via  employment  of  PDF  (Probability  Density  Function) 
approach  (see  Amin  and  Celik,  1999). 


tefi9(iQ 


Figure  8.14  Temperature  contours  at  15°  CA  atdc  for  3D  case  with  piston  bowl 


Figure  8.15  Temperature  contours  at  24°  CA  btdc  of  the  3D  case  with  2  valves 

At  the  second  stage  a  Cartesian  block-structured  grid  with  220,000  nodes  was  applied  to 
a  two-valve  flat-piston  geometry.  The  Smagorinsky  subgrid-scale  model  and  the  Magnussen 
eddy-dissipation  model  were  used  to  account  for  the  effects  of  turbulence  and  combustion 
respectively.  Figure  8.15  presents  the  temperature  contours  of  this  case  on  a  horizontal  cross- 
section  just  before  TDC.  Spray  distribution  and  dynamics  was  found  to  be  grid-sensitive.  A  finer 
grid  leads  usually  to  a  slower  spray  penetration  and  higher  breakup  frequencies.  Autoignition  of 
the  fuel  was  modeled  by  a  standard  Arrhenius  reaction  mechanism.  Time  step  was  varied 
depending  on  the  time-scale  of  the  limiting  process.  Thus,  in  the  evaporation  phase  (285°-360° 
CA)  the  time  steps  could  be  as  small  as  10'^  seconds  which  was  due  to  the  limits  imposed  by  the 
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evaporation  process  of  smallest  droplets.  In  this  case  the  geometry  is  essentially  asymmetric  with 
two  valves  located  off  the  symmetry  plane  of  the  cylinder,  but  the  jet  spray  was  axe-symmetric  at 
the  inlet  of  the  nozzle.  As  a  result  the  ignition  starts  asymmetrically  with  multiple  ignition 
regions  (Fig.8.15).  These  regions  are  confined  to  the  peripheral  part  of  the  cylinder.  Both 
turbulent  dynamics  and  the  geometry  may  contribute  to  the  symmetry  breaking  in  this  case. 


The  onset  of  combustion  in  both  cases  was  marked  by  an  increase  in  turbulent  intensity, 
which  was  more  pronounced  in  the  case  when  the  piston  bowl  was  present  (Fig.  8.16).  This 
result  is  supported  by  experimental  studies  by  Driscoll  et  al  (Driscoll  et  ah,  1986),  where  the 
turbulence  intensity  increase  associated  with  combustion  was  demonstrated  for  swirling  flows  of 
various  swirl  intensity.  This  finding  is  also  consistent  with  the  results  of  Rapid  Distortion  theory 
(Hoult  and  Wong,  1980).  A  more  visual  understanding  of  the  three-dimensional  nature  of 
combustion  the  iso-contours  of  temperature  during  combustion  are  presented  in  Figure  8.17. 
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b)  Two-valve  assembly 
Figure  8.16  Velocity  fluctuations  at  the  auto-ignition  point 


a)  Bowl  Geometry 


Figure  8.17  Temperature  Iso-Contours  during  Combustion  in  Bowl  Geometry 


105 


CONCLUSIONS 


A  systematic  investigation  of  numerical  prediction  of  in-cylinder  turbulence  (i.e. 
instantaneous  random  fluctuating  components  of  flow  variables)  in  internal  combustion  engines 
is  performed.  For  this  purpose,  the  well-known  KIVA  code  has  been  utilized.  The  results  of  the 
simulated  benchmark  cases  have  revealed  that,  though  KIVA  is  not  the  most  suitable  code  for 
LES,  nevertheless,  by  carefully  controlling  the  numerical  errors  and  using  relatively  fine  grid 
resolution  and  a  very  small  time  step  the  unsteady  dynamics  of  turbulent  flows  can  be  captured 
reasonably  well.  Particularly  the  time  averaged  predictions  for  the  free  swirling  annular  jet  are  in 
reasonable  agreement  with  measurements  and  with  results  obtained  from  another  well  tested  LES 
code,  which  gives  confidence  that  the  in-cylinder  turbulence  inside  engine  cylinders  can  be 
studied  with  this  code. 

The  present  engine  applications  show  that  the  instantaneous  large  scale  flow  structures 
(i.e.  coherent  large  eddies)  can  be  captured,  and  the  predicted  turbulence  statistics  are  in  good 
qualitative  agreement  with  the  experimental  data  in  a  similar  engine.  In  particular,  the  good 
correspondence  of  normalized  spectra  to  that  obtained  from  measurements  by  other  authors  is 
encouraging.  It  can  be  said  that  a  significant  portion  of  the  energy  spectra  can  be  resolved  within 
the  range  of  100  to  10,000  Hz.  However,  the  reader  should  be  cautioned  that  the  cycle-to-cycle 
variations,  which  are  usually  significant,  could  not  be  accounted  for  in  the  present  calculations. 
The  predicted  trends  of  velocity  fluctuations  during  the  intake  as  a  function  of  crank  angle  as 
well  as  its  magnitude  are  also  in  good  agreement  with  experimentally  observed  trends.  Close 
agreement  of  fluctuating  velocity  component  produced  on  different  grids  and  a  good  comparison 
with  experimental  data  indicate  that  the  numerical  diffusion  errors  can  not  be  a  dominant  factor 
in  the  defining  the  dynamics  of  the  resolved  flow-field. 

This  study  has  also  shown  that  contrary  to  some  arguments,  significant  turbulence  is 
likely  to  be  generated  by  a  carefully  designed  bowl,  and  that  the  intake  turbulence  decays  rapidly 
towards  the  TDC  during  compression  stroke.  Instability  seems  to  be  induced  by  the  unsteady 
flow  separation  at  the  edge  of  the  bowl  and  the  squish  region. 

It  has  been  further  demonstrated  that  the  statistics  of  coherent  large  flow  structures  can  be 
studied  using  the  widely  used  engine  code  KIVA  on  a  workstation  level.  This  opens  a  whole  new 
avenue  in  simulation  of  in-cylinder  turbulence  for  IC-engines.  Extension  of  this  study  to  cases 
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with  combustion  and  many  cycles  is  currently  underway.  The  preliminary  results  from  spray 
combustion  are  encouraging  in  this  regard. 

Simulations  of  the  combustion  stroke  show  that  the  intensity  of  velocity  fluctuations 
increases  at  the  onset  of  combustion,  though  at  a  level  considerably  lower  than  observed  during 
the  intake  stroke.  Temperature  distributions  captured  at  a  single  point  show  a  considerable 
intermittency  at  least  at  the  initial  stage  of  combustion. 

Because  of  the  general  decay  of  turbulence  after  the  intake  stroke,  shown  in  experiments 
and  confirmed  by  the  present  calculations,  turbulence  generated  at  the  onset  of  combustion  is  not 
greatly  affected  by  the  intake  turbulence.  This  may  provide  some  justification  for  setting  up 
approximate  initial  conditions  at  the  intake  stroke,  i.e.  without  modeling  a  complicated  valve 
assembly,  and  still  reproduce  the  main  dynamics  of  the  combustion  process.  Spray  evaporation 
time  constants  may  significantly  reduce  the  execution  speed  by  limiting  the  iteration  time-step. 
Spray  distribution  and  dynamics  show  some  grid-dependence.  An  improvement  of  the  spray 
model  would  be  desirable. 

The  case  studies  have  shown  that  the  intake  turbulence  generation  and  decay  can  be 
predicted  in  good  agreement  with  experiments.  Turbulence  generated  during  the  intake  decays 
rapidly,  and  the  turbulence  generated  by  the  bowl  is  relatively  small.  Turbulence  generation 
increases  with  the  onset  of  combustion. 

The  accuracy  and  efficiency  of  the  primary  simulation  code  KIVA-3  was  improved  in 
many  ways,  as  the  time  accuracy,  spatial  accuracy,  and  the  efficiency  of  the  pressure  solver  was 
improved.  The  overall  computational  performance  was  improved  significantly  with  the 
distributed-memory  implementation  of  KIVA-3  (KIVA-3/MPI).  Moreover,  a  beowulf  cluster  of 
20  alpha  processors  has  been  developed  for  complex  computations.  Twelve  nodes  are  made 
operational  with  the  MPI  (Message  Passing  Interface)  library,  and  successful  KIVA-3  runs  have 
been  performed  on  this  cluster.  Although  it  was  completed  rather  late  to  have  a  major  impact  on 
the  current  simulations  reported  on  here,  it  should  be  an  important  asset  for  future  simulations. 

Turbulence  measurement  capabilities  with  non-intrusive  Doppler  Global  velocimetry 
were  also  developed  which  should  have  an  impact  in  future  validation  studies. 
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RECOMMENDATIONS 


This  study  has  demonstrated  that  much  can  be  learned  from  large  eddy  simulations  of 
turbulent  flow  dynamics  and  combustion  in  diesel  combustion  chambers.  There  still  remain, 
however,  some  obstacles  to  be  overcome.  First  of  all  in  order  to  ensure  practicality  in 
engineering  applications,  the  turnaround  time  for  the  computations  should  be  decreased 
drastically.  This  can  be  accomplished  by  resorting  to  parallel  programming  with  the  appropriate 
numerical  algorithms  and  the  use  of  parallel  computers.  A  parallel  implementation  of  the  KJVA 
code  on  work  station  clusters  with  distributed  memory  seems  to  be  a  viable  alternative  and  cost 
effective  for  research  institutions  as  well  as  industry.  To  accomplish  this,  the  preliminary  KIVA- 
3/MPI  parallel  version  needs  to  be  extended  to  more  general  domains,  and  sub-models,  such  as 
combustion  and  spray  dynamics  need  to  be  included. 

Since  the  laminar  flame  thickness  is  very  thin  and  beyond  the  realm  of  LES  resolutions,  a 
robust,  well  tested,  and  validated  sub-grid  scale  combustion  model  is  essential  for  accurate 
prediction  of  diffusion  flames.  The  computations  of  Lagrangian  spray  dynamics  with  droplet 
evaporation  consume  most  of  the  CPU  time.  Again,  this  can  be  improved  by  implementing  a 
parallel  particle-tracking  algorithm,  but  it  is  also  crucial  to  improve  the  stability  limits  of  spray 
calculations  especially  when  the  spray  impinges  at  the  cylinder  walls. 

It  is  also  highly  desirable  to  obtain  a  good  set  of  experimental  data  specifically  for 
validation  of  large  eddy  simulations.  The  experiments  need  to  be  performed  in  conjunction  with 
the  computations,  and  cover  a  range  of  engine  speeds  with  and  without  combustion,  while 
keeping  the  valve  and  cylinder  geometry  relatively  simple. 
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